#include "PolynomialPack.h"

/* Common Block Declarations */

union {
    struct {
	doublereal break__[100], coef[2000];
	integer l, kpm;
    } _1;
    struct {
	doublereal break__[100], coef[2000];
	integer l, k;
    } _2;
} approx_;

#define approx_1 (approx_._1)
#define approx_2 (approx_._2)

struct {
    integer m, iside;
    doublereal xside[10];
} side_;

#define side_1 side_

struct {
    integer itermx, k;
    doublereal rho[19];
} other_;

#define other_1 other_

struct {
    integer ntau;
} i4data_;

#define i4data_1 i4data_

union {
    struct {
	doublereal tau[200], gtau[200], weight[200];
    } _1;
    struct {
	doublereal tau[200], gtau[200], weight[200], totalw;
    } _2;
} r8data_;

#define r8data_1 (r8data_._1)
#define r8data_2 (r8data_._2)

struct {
    doublereal size;
} rount_;

#define rount_1 rount_

/* Table of constant values */

static integer c__1 = 1;
static integer c__2 = 2;
static integer c__2000 = 2000;
static integer c__4 = 4;
static integer c__0 = 0;
static doublereal c_b213 = 1.;

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::banfac_(doublereal *w, integer *nroww, integer *nrow, 
	integer *nbandl, integer *nbandu, integer *iflag)
{
    /* System generated locals */
    integer w_dim1, w_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i__, j, k, ipk, jmax, kmax, midmk;
    static doublereal pivot;
    static integer nrowm1, middle;
    static doublereal factor;

/* *********************************************************************72 */

/* c BANFAC factors a banded matrix without pivoting. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  returns in  w  the lu-factorization (without pivoting) of the banded */
/*  matrix  a  of order  nrow  with  (nbandl + 1 + nbandu) bands or diag- */
/*  onals in the work array  w . */

/* ******  i n p u t  ****** */
/*  w.....work array of size  (nroww,nrow)  containing the interesting */
/*        part of a banded matrix  a , with the diagonals or bands of  a */
/*        stored in the rows of  w , while columns of  a  correspond to */
/*        columns of  w . this is the storage mode used in  linpack  and */
/*        results in efficient innermost loops. */
/*           explicitly,  a  has  nbandl  bands below the diagonal */
/*                            +     1     (main) diagonal */
/*                            +   nbandu  bands above the diagonal */
/*        and thus, with    middle = nbandu + 1, */
/*          a(i+j,j)  is in  w(i+middle,j)  for i=-nbandu,...,nbandl */
/*                                              j=1,...,nrow . */
/*        for example, the interesting entries of a (1,2)-banded matrix */
/*        of order  9  would appear in the first  1+1+2 = 4  rows of  w */
/*        as follows. */





/*        all other entries of  w  not identified in this way with an en- */
/*        try of  a  are never referenced . */
/*  nroww.....row dimension of the work array  w . */
/*        must be  .ge.  nbandl + 1 + nbandu  . */
/*  nbandl.....number of bands of  a  below the main diagonal */
/*  nbandu.....number of bands of  a  above the main diagonal . */

/* ******  o u t p u t  ****** */
/*  iflag.....integer indicating success( = 1) or failure ( = 2) . */
/*     if  iflag = 1, then */
/*  w.....contains the lu-factorization of  a  into a unit lower triangu- */
/*        lar matrix  l  and an upper triangular matrix  u (both banded) */
/*        and stored in customary fashion over the corresponding entries */
/*        of  a . this makes it possible to solve any particular linear */
/*        system  a*x = b  for  x  by a */
/*              call banslv ( w, nroww, nrow, nbandl, nbandu, b ) */
/*        with the solution x  contained in  b  on return . */
/*     if  iflag = 2, then */
/*        one of  nrow-1, nbandl,nbandu failed to be nonnegative, or else */
/*        one of the potential pivots was found to be zero indicating */
/*        that  a  does not have an lu-factorization. this implies that */
/*        a  is singular in case it is totally positive . */

/* ******  m e t h o d  ****** */
/*     gauss elimination  w i t h o u t  pivoting is used. the routine is */
/*  intended for use with matrices  a  which do not require row inter- */
/*  changes during factorization, especially for the  t o t a l l y */
/*  p o s i t i v e  matrices which occur in spline calculations. */
/*     the routine should not be used for an arbitrary banded matrix. */


    /* Parameter adjustments */
    w_dim1 = *nroww;
    w_offset = 1 + w_dim1;
    w -= w_offset;

    /* Function Body */
    *iflag = 1;
    middle = *nbandu + 1;
/*                         w(middle,.) contains the main diagonal of  a . */
    nrowm1 = *nrow - 1;
    if (nrowm1 < 0) {
	goto L999;
    } else if (nrowm1 == 0) {
	goto L900;
    } else {
	goto L1;
    }
L1:
    if (*nbandl > 0) {
	goto L10;
    }
/*                a is upper triangular. check that diagonal is nonzero . */
    i__1 = nrowm1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	if (w[middle + i__ * w_dim1] == 0.) {
	    goto L999;
	}
/* L5: */
    }
    goto L900;
L10:
    if (*nbandu > 0) {
	goto L20;
    }
/*              a is lower triangular. check that diagonal is nonzero and */
/*                 divide each column by its diagonal . */
    i__1 = nrowm1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	pivot = w[middle + i__ * w_dim1];
	if (pivot == 0.) {
	    goto L999;
	}
/* Computing MIN */
	i__2 = *nbandl, i__3 = *nrow - i__;
	jmax = min(i__2,i__3);
	i__2 = jmax;
	for (j = 1; j <= i__2; ++j) {
/* L15: */
	    w[middle + j + i__ * w_dim1] /= pivot;
	}
    }
    return 0;

/*        a  is not just a triangular matrix. construct lu factorization */
L20:
    i__2 = nrowm1;
    for (i__ = 1; i__ <= i__2; ++i__) {
/*                                  w(middle,i)  is pivot for i-th step . */
	pivot = w[middle + i__ * w_dim1];
	if (pivot == 0.) {
	    goto L999;
	}
/*                 jmax  is the number of (nonzero) entries in column  i */
/*                     below the diagonal . */
/* Computing MIN */
	i__1 = *nbandl, i__3 = *nrow - i__;
	jmax = min(i__1,i__3);
/*              divide each entry in column  i  below diagonal by pivot . */
	i__1 = jmax;
	for (j = 1; j <= i__1; ++j) {
/* L32: */
	    w[middle + j + i__ * w_dim1] /= pivot;
	}
/*                 kmax  is the number of (nonzero) entries in row  i  to */
/*                     the right of the diagonal . */
/* Computing MIN */
	i__1 = *nbandu, i__3 = *nrow - i__;
	kmax = min(i__1,i__3);
/*                  subtract  a(i,i+k)*(i-th column) from (i+k)-th column */
/*                  (below row  i ) . */
	i__1 = kmax;
	for (k = 1; k <= i__1; ++k) {
	    ipk = i__ + k;
	    midmk = middle - k;
	    factor = w[midmk + ipk * w_dim1];
	    i__3 = jmax;
	    for (j = 1; j <= i__3; ++j) {
/* L40: */
		w[midmk + j + ipk * w_dim1] -= w[middle + j + i__ * w_dim1] * 
			factor;
	    }
	}
/* L50: */
    }
/*                                       check the last diagonal entry . */
L900:
    if (w[middle + *nrow * w_dim1] != 0.) {
	return 0;
    }
L999:
    *iflag = 2;
    return 0;
} /* banfac_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::banslv_(doublereal *w, integer *nroww, integer *nrow, 
	integer *nbandl, integer *nbandu, doublereal *b)
{
    /* System generated locals */
    integer w_dim1, w_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i__, j, jmax, nrowm1, middle;

/* *********************************************************************72 */

/* c BANSLV solves a banded linear system A * X = B factored by BANFAC. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  companion routine to  banfac . it returns the solution  x  of the */
/*  linear system  a*x = b  in place of  b , given the lu-factorization */
/*  for  a  in the workarray  w . */

/* ******  i n p u t  ****** */
/*  w, nroww,nrow,nbandl,nbandu.....describe the lu-factorization of a */
/*        banded matrix  a  of order  nrow  as constructed in  banfac . */
/*        for details, see  banfac . */
/*  b.....right side of the system to be solved . */

/* ******  o u t p u t  ****** */
/*  b.....contains the solution  x , of order  nrow . */

/* ******  m e t h o d  ****** */
/*     (with  a = l*u, as stored in  w,) the unit lower triangular system */
/*  l(u*x) = b  is solved for  y = u*x, and  y  stored in  b . then the */
/*  upper triangular system  u*x = y  is solved for  x  . the calcul- */
/*  ations are so arranged that the innermost loops stay within columns. */

    /* Parameter adjustments */
    --b;
    w_dim1 = *nroww;
    w_offset = 1 + w_dim1;
    w -= w_offset;

    /* Function Body */
    middle = *nbandu + 1;
    if (*nrow == 1) {
	goto L49;
    }
    nrowm1 = *nrow - 1;
    if (*nbandl == 0) {
	goto L30;
    }
/*                                 forward pass */
/*            for i=1,2,...,nrow-1, subtract  right side(i)*(i-th column */
/*            of  l )  from right side  (below i-th row) . */
    i__1 = nrowm1;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing MIN */
	i__2 = *nbandl, i__3 = *nrow - i__;
	jmax = min(i__2,i__3);
	i__2 = jmax;
	for (j = 1; j <= i__2; ++j) {
/* L21: */
	    b[i__ + j] -= b[i__] * w[middle + j + i__ * w_dim1];
	}
    }
/*                                 backward pass */
/*            for i=nrow,nrow-1,...,1, divide right side(i) by i-th diag- */
/*            onal entry of  u, then subtract  right side(i)*(i-th column */
/*            of  u)  from right side  (above i-th row). */
L30:
    if (*nbandu > 0) {
	goto L40;
    }
/*                                a  is lower triangular . */
    i__2 = *nrow;
    for (i__ = 1; i__ <= i__2; ++i__) {
/* L31: */
	b[i__] /= w[i__ * w_dim1 + 1];
    }
    return 0;
L40:
    i__ = *nrow;
L41:
    b[i__] /= w[middle + i__ * w_dim1];
/* Computing MIN */
    i__2 = *nbandu, i__1 = i__ - 1;
    jmax = min(i__2,i__1);
    i__2 = jmax;
    for (j = 1; j <= i__2; ++j) {
/* L45: */
	b[i__ - j] -= b[i__] * w[middle - j + i__ * w_dim1];
    }
    --i__;
    if (i__ > 1) {
	goto L41;
    }
L49:
    b[1] /= w[middle + w_dim1];
    return 0;
} /* banslv_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::bchfac_(doublereal *w, integer *nbands, integer *nrow, 
	doublereal *diag)
{
    /* System generated locals */
    integer w_dim1, w_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i__, j, n, imax, jmax;
    static doublereal ratio;

/* *********************************************************************72 */

/* c BCHFAC constructs a Cholesky factorization of a matrix. */

/*  from  * a practical guide to splines *  by c. de boor */
/* onstructs cholesky factorization */
/*                     c  =  l * d * l-transpose */
/*  with l unit lower triangular and d diagonal, for given matrix c of */
/*  order  n r o w , in case  c  is (symmetric) positive semidefinite */
/*  and  b a n d e d , having  n b a n d s  diagonals at and below the */
/*  main diagonal. */

/* ******  i n p u t  ****** */
/*  nrow.....is the order of the matrix  c . */
/*  nbands.....indicates its bandwidth, i.e., */
/*          c(i,j) = 0 for i-j .ge. nbands . */
/*  w.....workarray of size (nbands,nrow)  containing the  nbands  diago- */
/*        nals in its rows, with the main diagonal in row  1 . precisely, */
/*        w(i,j)  contains  c(i+j-1,j), i=1,...,nbands, j=1,...,nrow. */
/*          for example, the interesting entries of a seven diagonal sym- */
/*        metric matrix  c  of order  9  would be stored in  w  as */






/*        all other entries of  w  not identified in this way with an en- */
/*        try of  c  are never referenced . */
/*  diag.....is a work array of length  nrow . */

/* ******  o u t p u t  ****** */
/*  w.....contains the cholesky factorization  c = l*d*l-transp, with */
/*        w(1,i) containing  1/d(i,i) */
/*        and  w(i,j)  containing  l(i-1+j,j), i=2,...,nbands. */

/* ******  m e t h o d  ****** */
/*   gauss elimination, adapted to the symmetry and bandedness of  c , is */
/*   used . */
/*     near zero pivots are handled in a special way. the diagonal ele- */
/*  ment c(n,n) = w(1,n) is saved initially in  diag(n), all n. at the n- */
/*  th elimination step, the current pivot element, viz.  w(1,n), is com- */
/*  pared with its original value, diag(n). if, as the result of prior */
/*  elimination steps, this element has been reduced by about a word */
/*  length, (i.e., if w(1,n)+diag(n) .le. diag(n)), then the pivot is de- */
/*  clared to be zero, and the entire n-th row is declared to be linearly */
/*  dependent on the preceding rows. this has the effect of producing */
/*   x(n) = 0  when solving  c*x = b  for  x, regardless of  b. justific- */
/*  ation for this is as follows. in contemplated applications of this */
/*  program, the given equations are the normal equations for some least- */
/*  squares approximation problem, diag(n) = c(n,n) gives the norm-square */
/*  of the n-th basis function, and, at this point,  w(1,n)  contains the */
/*  norm-square of the error in the least-squares approximation to the n- */
/*  th basis function by linear combinations of the first n-1 . having */
/*  w(1,n)+diag(n) .le. diag(n) signifies that the n-th function is lin- */
/*  early dependent to machine accuracy on the first n-1 functions, there */
/*  fore can safely be left out from the basis of approximating functions */
/*     the solution of a linear system */
/*                       c*x = b */
/*   is effected by the succession of the following  t w o  calls: */
/*     call bchfac ( w, nbands, nrow, diag )       , to get factorization */
/*     call bchslv ( w, nbands, nrow, b )          , to solve for x. */

    /* Parameter adjustments */
    --diag;
    w_dim1 = *nbands;
    w_offset = 1 + w_dim1;
    w -= w_offset;

    /* Function Body */
    if (*nrow > 1) {
	goto L9;
    }
    if (w[w_dim1 + 1] > 0.) {
	w[w_dim1 + 1] = 1. / w[w_dim1 + 1];
    }
    return 0;
/*                                        store diagonal of  c  in  diag. */
L9:
    i__1 = *nrow;
    for (n = 1; n <= i__1; ++n) {
/* L10: */
	diag[n] = w[n * w_dim1 + 1];
    }
/*                                                        factorization . */
    i__1 = *nrow;
    for (n = 1; n <= i__1; ++n) {
	if (w[n * w_dim1 + 1] + diag[n] > diag[n]) {
	    goto L15;
	}
	i__2 = *nbands;
	for (j = 1; j <= i__2; ++j) {
/* L14: */
	    w[j + n * w_dim1] = 0.;
	}
	goto L20;
L15:
	w[n * w_dim1 + 1] = (float)1. / w[n * w_dim1 + 1];
/* Computing MIN */
	i__2 = *nbands - 1, i__3 = *nrow - n;
	imax = min(i__2,i__3);
	if (imax < 1) {
	    goto L20;
	}
	jmax = imax;
	i__2 = imax;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    ratio = w[i__ + 1 + n * w_dim1] * w[n * w_dim1 + 1];
	    i__3 = jmax;
	    for (j = 1; j <= i__3; ++j) {
/* L17: */
		w[j + (n + i__) * w_dim1] -= w[j + i__ + n * w_dim1] * ratio;
	    }
	    --jmax;
/* L18: */
	    w[i__ + 1 + n * w_dim1] = ratio;
	}
L20:
	;
    }
    return 0;
} /* bchfac_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::bchslv_(doublereal *w, integer *nbands, integer *nrow, 
	doublereal *b)
{
    /* System generated locals */
    integer w_dim1, w_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer j, n, jmax, nbndm1;

/* *********************************************************************72 */

/* c BCHSLV solves a banded symmetric positive definite system. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  solves the linear system     c*x = b   of order  n r o w  for  x */
/*  provided  w  contains the cholesky factorization for the banded (sym- */
/*  metric) positive definite matrix  c  as constructed in the subroutine */
/*    b c h f a c  (quo vide). */

/* ******  i n p u t  ****** */
/*  nrow.....is the order of the matrix  c . */
/*  nbands.....indicates the bandwidth of  c . */
/*  w.....contains the cholesky factorization for  c , as output from */
/*        subroutine bchfac  (quo vide). */
/*  b.....the vector of length  n r o w  containing the right side. */

/* ******  o u t p u t  ****** */
/*  b.....the vector of length  n r o w  containing the solution. */

/* ******  m e t h o d  ****** */
/*  with the factorization  c = l*d*l-transpose  available, where  l  is */
/*  unit lower triangular and  d  is diagonal, the triangular system */
/*  l*y = b  is solved for  y (forward substitution), y is stored in  b, */
/*  the vector  d**(-1)*y is computed and stored in  b, then the triang- */
/*  ular system  l-transpose*x = d**(-1)*y is solved for  x (backsubstit- */
/*  ution). */

    /* Parameter adjustments */
    --b;
    w_dim1 = *nbands;
    w_offset = 1 + w_dim1;
    w -= w_offset;

    /* Function Body */
    if (*nrow > 1) {
	goto L21;
    }
    b[1] *= w[w_dim1 + 1];
    return 0;

/*                 forward substitution. solve l*y = b for y, store in b. */
L21:
    nbndm1 = *nbands - 1;
    i__1 = *nrow;
    for (n = 1; n <= i__1; ++n) {
/* Computing MIN */
	i__2 = nbndm1, i__3 = *nrow - n;
	jmax = min(i__2,i__3);
	if (jmax < 1) {
	    goto L30;
	}
	i__2 = jmax;
	for (j = 1; j <= i__2; ++j) {
/* L25: */
	    b[j + n] -= w[j + 1 + n * w_dim1] * b[n];
	}
L30:
	;
    }

/*     backsubstitution. solve l-transp.x = d**(-1)*y  for x, store in b. */
    n = *nrow;
L39:
    b[n] *= w[n * w_dim1 + 1];
/* Computing MIN */
    i__1 = nbndm1, i__2 = *nrow - n;
    jmax = min(i__1,i__2);
    if (jmax < 1) {
	goto L40;
    }
    i__1 = jmax;
    for (j = 1; j <= i__1; ++j) {
/* L35: */
	b[n] -= w[j + 1 + n * w_dim1] * b[j + n];
    }
L40:
    --n;
    if (n > 0) {
	goto L39;
    }
    return 0;
} /* bchslv_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::bsplpp_(doublereal *t, doublereal *bcoef, integer *n, 
	integer *k, doublereal *scrtch, doublereal *break__, doublereal *coef,
	 integer *l)
{
    /* System generated locals */
    integer coef_dim1, coef_offset, scrtch_dim1, scrtch_offset, i__1, i__2, 
	    i__3;

    /* Local variables */
    static integer i__, j, jp1, kmj;
    static doublereal sum, diff;
    static integer left;
    static doublereal biatx[20], factor;
    static integer lsofar;

/* *********************************************************************72 */

/* c BSPLPP converts from B-spline to piecewise polynomial form. */

/*  from  * a practical guide to splines *  by c. de Boor (7 may 92) */
/* alls  bsplvb */

/* onverts the b-representation  t, bcoef, n, k  of some spline into its */
/*  pp-representation  break, coef, l, k . */

/* ******  i n p u t  ****** */
/*  t.....knot sequence, of length  n+k */
/*  bcoef.....b-spline coefficient sequence, of length  n */
/*  n.....length of  bcoef  and  dimension of spline space  spline(k,t) */
/*  k.....order of the spline */

/*  w a r n i n g  . . .  the restriction   k .le. kmax (= 20)   is impo- */
/*        sed by the arbitrary dimension statement for  biatx  below, but */
/*        is  n o w h e r e   c h e c k e d   for. */

/* ******  w o r k   a r e a  ****** */
/*  scrtch......of size  (k,k) , needed to contain bcoeffs of a piece of */
/*        the spline and its  k-1  derivatives */

/* ******  o u t p u t  ****** */
/*  break.....breakpoint sequence, of length  l+1, contains (in increas- */
/*        ing order) the distinct points in the sequence  t(k),...,t(n+1) */
/*  coef.....array of size (k,l), with  coef(i,j) = (i-1)st derivative of */
/*        spline at break(j) from the right */
/*  l.....number of polynomial pieces which make up the spline in the in- */
/*        terval  (t(k), t(n+1)) */

/* ******  m e t h o d  ****** */
/*     for each breakpoint interval, the  k  relevant b-coeffs of the */
/*  spline are found and then differenced repeatedly to get the b-coeffs */
/*  of all the derivatives of the spline on that interval. the spline and */
/*  its first  k-1  derivatives are then evaluated at the left end point */
/*  of that interval, using  bsplvb  repeatedly to obtain the values of */
/*  all b-splines of the appropriate order at that point. */


    /* Parameter adjustments */
    --bcoef;
    scrtch_dim1 = *k;
    scrtch_offset = 1 + scrtch_dim1;
    scrtch -= scrtch_offset;
    --t;
    coef_dim1 = *k;
    coef_offset = 1 + coef_dim1;
    coef -= coef_offset;
    --break__;

    /* Function Body */
    lsofar = 0;
    break__[1] = t[*k];
    i__1 = *n;
    for (left = *k; left <= i__1; ++left) {
/*                                find the next nontrivial knot interval. */
	if (t[left + 1] == t[left]) {
	    goto L50;
	}
	++lsofar;
	break__[lsofar + 1] = t[left + 1];
	if (*k > 1) {
	    goto L9;
	}
	coef[lsofar * coef_dim1 + 1] = bcoef[left];
	goto L50;
/*        store the k b-spline coeff.s relevant to current knot interval */
/*                             in  scrtch(.,1) . */
L9:
	i__2 = *k;
	for (i__ = 1; i__ <= i__2; ++i__) {
/* L10: */
	    scrtch[i__ + scrtch_dim1] = bcoef[left - *k + i__];
	}

/*        for j=1,...,k-1, compute the  k-j  b-spline coeff.s relevant to */
/*        current knot interval for the j-th derivative by differencing */
/*        those for the (j-1)st derivative, and store in scrtch(.,j+1) . */
	i__2 = *k;
	for (jp1 = 2; jp1 <= i__2; ++jp1) {
	    j = jp1 - 1;
	    kmj = *k - j;
	    i__3 = kmj;
	    for (i__ = 1; i__ <= i__3; ++i__) {
		diff = t[left + i__] - t[left + i__ - kmj];
		if (diff > 0.) {
		    scrtch[i__ + jp1 * scrtch_dim1] = (scrtch[i__ + 1 + j * 
			    scrtch_dim1] - scrtch[i__ + j * scrtch_dim1]) / 
			    diff;
		}
/* L20: */
	    }
	}

/*        for  j = 0, ..., k-1, find the values at  t(left)  of the  j+1 */
/*        b-splines of order  j+1  whose support contains the current */
/*        knot interval from those of order  j  (in  biatx ), then comb- */
/*        ine with the b-spline coeff.s (in scrtch(.,k-j) ) found earlier */
/*        to compute the (k-j-1)st derivative at  t(left)  of the given */
/*        spline. */
/*           note. if the repeated calls to  bsplvb  are thought to gene- */
/*        rate too much overhead, then replace the first call by */
/*           biatx(1) = 1. */
/*        and the subsequent call by the statement */
/*           j = jp1 - 1 */
/*        followed by a direct copy of the lines */
/*           deltar(j) = t(left+j) - x */
/*                  ...... */
/*           biatx(j+1) = saved */
/*        from  bsplvb . deltal(kmax)  and  deltar(kmax)  would have to */
/*        appear in a dimension statement, of course. */

	bsplvb_(&t[1], &c__1, &c__1, &t[left], &left, biatx);
	coef[*k + lsofar * coef_dim1] = scrtch[*k * scrtch_dim1 + 1];
	i__3 = *k;
	for (jp1 = 2; jp1 <= i__3; ++jp1) {
	    bsplvb_(&t[1], &jp1, &c__2, &t[left], &left, biatx);
	    kmj = *k + 1 - jp1;
	    sum = 0.;
	    i__2 = jp1;
	    for (i__ = 1; i__ <= i__2; ++i__) {
/* L28: */
		sum = biatx[i__ - 1] * scrtch[i__ + kmj * scrtch_dim1] + sum;
	    }
/* L30: */
	    coef[kmj + lsofar * coef_dim1] = sum;
	}
L50:
	;
    }
    *l = lsofar;
    if (*k == 1) {
	return 0;
    }
    factor = 1.;
    i__1 = *k;
    for (i__ = 2; i__ <= i__1; ++i__) {
	factor *= (doublereal) (*k + 1 - i__);
	i__3 = lsofar;
	for (j = 1; j <= i__3; ++j) {
/* L60: */
	    coef[i__ + j * coef_dim1] *= factor;
	}
    }
    return 0;
} /* bsplpp_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::bsplvb_(doublereal *t, integer *jhigh, integer *index, 
	doublereal *x, integer *left, doublereal *biatx)
{
    /* Initialized data */

    static integer j = 1;

    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i__, jp1;
    static doublereal term, saved, deltal[20], deltar[20];

/* *********************************************************************72 */

/* c BSPLVB evaluates B-splines at a point X with a given knot sequence. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alculates the value of all possibly nonzero b-splines at  x  of order */

/*               jout  =  max( jhigh , (j+1)*(index-1) ) */

/*  with knot sequence  t . */

/* ******  i n p u t  ****** */
/*  t.....knot sequence, of length  left + jout  , assumed to be nonde- */
/*        creasing.  a s s u m p t i o n . . . . */
/*                       t(left)  .lt.  t(left + 1)   . */
/*   d i v i s i o n  b y  z e r o  will result if  t(left) = t(left+1) */
/*  jhigh, */
/*  index.....integers which determine the order  jout = max(jhigh, */
/*        (j+1)*(index-1))  of the b-splines whose values at  x  are to */
/*        be returned.  index  is used to avoid recalculations when seve- */
/*        ral columns of the triangular array of b-spline values are nee- */
/*        ded (e.g., in  bsplpp  or in  bsplvd ). precisely, */
/*                     if  index = 1 , */
/*        the calculation starts from scratch and the entire triangular */
/*        array of b-spline values of orders 1,2,...,jhigh  is generated */
/*        order by order , i.e., column by column . */
/*                     if  index = 2 , */
/*        only the b-spline values of order  j+1, j+2, ..., jout  are ge- */
/*        nerated, the assumption being that  biatx , j , deltal , deltar */
/*        are, on entry, as they were on exit at the previous call. */
/*           in particular, if  jhigh = 0, then  jout = j+1, i.e., just */
/*        the next column of b-spline values is generated. */

/*  w a r n i n g . . .  the restriction   jout .le. jmax (= 20)  is im- */
/*        posed arbitrarily by the dimension statement for  deltal  and */
/*        deltar  below, but is  n o w h e r e  c h e c k e d  for . */

/*  x.....the point at which the b-splines are to be evaluated. */
/*  left.....an integer chosen (usually) so that */
/*                  t(left) .le. x .le. t(left+1)  . */

/* ******  o u t p u t  ****** */
/*  biatx.....array of length  jout , with  biatx(i)  containing the val- */
/*        ue at  x  of the polynomial of order  jout  which agrees with */
/*        the b-spline  b(left-jout+i,jout,t)  on the interval (t(left), */
/*        t(left+1)) . */

/* ******  m e t h o d  ****** */
/*  the recurrence relation */

/*                       x - t(i)              t(i+j+1) - x */
/*     b(i,j+1)(x)  =  -----------b(i,j)(x) + ---------------b(i+1,j)(x) */
/*                     t(i+j)-t(i)            t(i+j+1)-t(i+1) */

/*  is used (repeatedly) to generate the (j+1)-vector  b(left-j,j+1)(x), */
/*  ...,b(left,j+1)(x)  from the j-vector  b(left-j+1,j)(x),..., */
/*  b(left,j)(x), storing the new values in  biatx  over the old. the */
/*  facts that */
/*            b(i,1) = 1  if  t(i) .le. x .lt. t(i+1) */
/*  and that */
/*            b(i,j)(x) = 0  unless  t(i) .le. x .lt. t(i+j) */
/*  are used. the particular organization of the calculations follows al- */
/*  gorithm  (8)  in chapter x of the text. */


/*     dimension biatx(jout), t(left+jout) */
/* urrent fortran standard makes it impossible to specify the length of */
/*  t  and of  biatx  precisely without the introduction of otherwise */
/*  superfluous additional arguments. */
    /* Parameter adjustments */
    --t;
    --biatx;

    /* Function Body */

    switch (*index) {
	case 1:  goto L10;
	case 2:  goto L20;
    }
L10:
    j = 1;
    biatx[1] = 1.;
    if (j >= *jhigh) {
	goto L99;
    }

L20:
    jp1 = j + 1;
    deltar[j - 1] = t[*left + j] - *x;
    deltal[j - 1] = *x - t[*left + 1 - j];
    saved = 0.;
    i__1 = j;
    for (i__ = 1; i__ <= i__1; ++i__) {
	term = biatx[i__] / (deltar[i__ - 1] + deltal[jp1 - i__ - 1]);
	biatx[i__] = saved + deltar[i__ - 1] * term;
/* L26: */
	saved = deltal[jp1 - i__ - 1] * term;
    }
    biatx[jp1] = saved;
    j = jp1;
    if (j < *jhigh) {
	goto L20;
    }

L99:
    return 0;
} /* bsplvb_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::bsplvd_(doublereal *t, integer *k, doublereal *x, 
	integer *left, doublereal *a, doublereal *dbiatx, integer *nderiv)
{
    /* System generated locals */
    integer a_dim1, a_offset, dbiatx_dim1, dbiatx_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i__, j, m, il, kp1;
    static doublereal sum;
    static integer jlow, kp1mm, mhigh, jp1mid;
    static doublereal fkp1mm, factor;
    static integer ideriv;
    static integer ldummy;

/* *********************************************************************72 */

/* c BSPLVD calculates the nonvanishing B-splines and derivatives at X. */

/*  from  * a practical guide to splines *  by c. de Boor (7 may 92) */
/* alls bsplvb */
/* alculates value and deriv.s of all b-splines which do not vanish at x */

/* ******  i n p u t  ****** */
/*  t     the knot array, of length left+k (at least) */
/*  k     the order of the b-splines to be evaluated */
/*  x     the point at which these values are sought */
/*  left  an integer indicating the left endpoint of the interval of */
/*        interest. the  k  b-splines whose support contains the interval */
/*               (t(left), t(left+1)) */
/*        are to be considered. */
/*  a s s u m p t i o n  - - -  it is assumed that */
/*               t(left) .lt. t(left+1) */
/*        division by zero will result otherwise (in  b s p l v b ). */
/*        also, the output is as advertised only if */
/*               t(left) .le. x .le. t(left+1) . */
/*  nderiv   an integer indicating that values of b-splines and their */
/*        derivatives up to but not including the  nderiv-th  are asked */
/*        for. ( nderiv  is replaced internally by the integer  m h i g h */
/*        in  (1,k)  closest to it.) */

/* ******  w o r k   a r e a  ****** */
/*  a     an array of order (k,k), to contain b-coeff.s of the derivat- */
/*        ives of a certain order of the  k  b-splines of interest. */

/* ******  o u t p u t  ****** */
/*  dbiatx   an array of order (k,nderiv). its entry  (i,m)  contains */
/*        value of  (m-1)st  derivative of  (left-k+i)-th  b-spline of */
/*        order  k  for knot sequence  t , i=1,...,k, m=1,...,nderiv. */

/* ******  m e t h o d  ****** */
/*  values at  x  of all the relevant b-splines of order k,k-1,..., */
/*  k+1-nderiv  are generated via  bsplvb  and stored temporarily in */
/*  dbiatx .  then, the b-coeffs of the required derivatives of the b- */
/*  splines of interest are generated by differencing, each from the pre- */
/*  ceding one of lower order, and combined with the values of b-splines */
/*  of corresponding order in  dbiatx  to produce the desired values . */

    /* Parameter adjustments */
    --t;
    a_dim1 = *k;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    dbiatx_dim1 = *k;
    dbiatx_offset = 1 + dbiatx_dim1;
    dbiatx -= dbiatx_offset;

    /* Function Body */
/* Computing MAX */
    i__1 = min(*nderiv,*k);
    mhigh = max(i__1,1);
/*     mhigh is usually equal to nderiv. */
    kp1 = *k + 1;
    i__1 = kp1 - mhigh;
    bsplvb_(&t[1], &i__1, &c__1, x, left, &dbiatx[dbiatx_offset]);
    if (mhigh == 1) {
	goto L99;
    }
/*     the first column of  dbiatx  always contains the b-spline values */
/*     for the current order. these are stored in column k+1-current */
/*     order  before  bsplvb  is called to put values for the next */
/*     higher order on top of it. */
    ideriv = mhigh;
    i__1 = mhigh;
    for (m = 2; m <= i__1; ++m) {
	jp1mid = 1;
	i__2 = *k;
	for (j = ideriv; j <= i__2; ++j) {
	    dbiatx[j + ideriv * dbiatx_dim1] = dbiatx[jp1mid + dbiatx_dim1];
/* L11: */
	    ++jp1mid;
	}
	--ideriv;
	i__2 = kp1 - ideriv;
	bsplvb_(&t[1], &i__2, &c__2, x, left, &dbiatx[dbiatx_offset]);
/* L15: */
    }

/*     at this point,  b(left-k+i, k+1-j)(x) is in  dbiatx(i,j) for */
/*     i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the */
/*     first column of  dbiatx  is already in final form. to obtain cor- */
/*     responding derivatives of b-splines in subsequent columns, gene- */
/*     rate their b-repr. by differencing, then evaluate at  x. */

    jlow = 1;
    i__1 = *k;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = *k;
	for (j = jlow; j <= i__2; ++j) {
/* L19: */
	    a[j + i__ * a_dim1] = 0.;
	}
	jlow = i__;
/* L20: */
	a[i__ + i__ * a_dim1] = 1.;
    }
/*     at this point, a(.,j) contains the b-coeffs for the j-th of the */
/*     k  b-splines of interest here. */

    i__1 = mhigh;
    for (m = 2; m <= i__1; ++m) {
	kp1mm = kp1 - m;
	fkp1mm = (doublereal) kp1mm;
	il = *left;
	i__ = *k;

/*        for j=1,...,k, construct b-coeffs of  (m-1)st  derivative of */
/*        b-splines from those for preceding derivative by differencing */
/*        and store again in  a(.,j) . the fact that  a(i,j) = 0  for */
/*        i .lt. j  is used. */
	i__2 = kp1mm;
	for (ldummy = 1; ldummy <= i__2; ++ldummy) {
	    factor = fkp1mm / (t[il + kp1mm] - t[il]);
/*           the assumption that t(left).lt.t(left+1) makes denominator */
/*           in  factor  nonzero. */
	    i__3 = i__;
	    for (j = 1; j <= i__3; ++j) {
/* L24: */
		a[i__ + j * a_dim1] = (a[i__ + j * a_dim1] - a[i__ - 1 + j * 
			a_dim1]) * factor;
	    }
	    --il;
/* L25: */
	    --i__;
	}

/*        for i=1,...,k, combine b-coeffs a(.,i) with b-spline values */
/*        stored in dbiatx(.,m) to get value of  (m-1)st  derivative of */
/*        i-th b-spline (of interest here) at  x , and store in */
/*        dbiatx(i,m). storage of this value over the value of a b-spline */
/*        of order m there is safe since the remaining b-spline derivat- */
/*        ives of the same order do not use this value due to the fact */
/*        that  a(j,i) = 0  for j .lt. i . */
	i__2 = *k;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    sum = 0.;
	    jlow = max(i__,m);
	    i__3 = *k;
	    for (j = jlow; j <= i__3; ++j) {
/* L35: */
		sum = a[j + i__ * a_dim1] * dbiatx[j + m * dbiatx_dim1] + sum;
	    }
/* L40: */
	    dbiatx[i__ + m * dbiatx_dim1] = sum;
	}
    }
L99:
    return 0;
} /* bsplvd_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::bspp2d_(doublereal *t, doublereal *bcoef, integer *n, 
	integer *k, integer *m, doublereal *scrtch, doublereal *break__, 
	doublereal *coef)
{
    /* System generated locals */
    integer bcoef_dim1, bcoef_offset, coef_dim1, coef_dim2, coef_offset, 
	    scrtch_dim1, scrtch_dim2, scrtch_offset, i__1, i__2, i__3, i__4;

    /* Local variables */
    static integer i__, j, mm, jp1, kmj;
    static doublereal sum, diff, fkmj;
    static integer left;
    static doublereal biatx[20];
    static integer lsofar;

/* *********************************************************************72 */

/* c BSPP2D converts from B-spline to piecewise polynomial representation. */

/*  from  * a practical guide to splines *  by c. de Boor */
/* alls  bsplvb */
/*  this is an extended version of  bsplpp  for use with tensor products -------- */

/* onverts the b-representation  t, bcoef(.,j), n, k  of some spline into -------- */
/*  its pp-representation  break, coef(j,.,.), l, k , j=1, ..., m  .     -------- */

/* ******  i n p u t  ****** */
/*  t.....knot sequence, of length  n+k */
/*  bcoef(.,j) b-spline coefficient sequence, of length  n ,j=1,...,m    -------- */
/*  n.....length of  bcoef  and  dimension of spline space  spline(k,t) */
/*  k.....order of the spline */

/*  w a r n i n g  . . .  the restriction   k .le. kmax (= 20)   is impo- */
/*        sed by the arbitrary dimension statement for  biatx  below, but */
/*        is  n o w h e r e   c h e c k e d   for. */

/*  m     number of data sets                                            ******** */

/* ******  w o r k   a r e a  ****** */
/*  scrtch   of size  (k,k,m), needed to contain bcoeffs of a piece of   ******** */
/*        the spline and its  k-1  derivatives   for each of the m sets  -------- */

/* ******  o u t p u t  ****** */
/*  break.....breakpoint sequence, of length  l+1, contains (in increas- */
/*        ing order) the distinct points in the sequence  t(k),...,t(n+1) */
/*  coef(mm,.,.)  array of size (k,n), with  coef(mm,i,j) = (i-1)st der- ******** */
/*        ivative of  mm-th  spline at break(j) from the right, mm=1,.,m ******** */

/* ******  m e t h o d  ****** */
/*     for each breakpoint interval, the  k  relevant b-coeffs of the */
/*  spline are found and then differenced repeatedly to get the b-coeffs */
/*  of all the derivatives of the spline on that interval. the spline and */
/*  its first  k-1  derivatives are then evaluated at the left end point */
/*  of that interval, using  bsplvb  repeatedly to obtain the values of */
/*  all b-splines of the appropriate order at that point. */


    /* Parameter adjustments */
    --break__;
    --t;
    coef_dim1 = *m;
    coef_dim2 = *k;
    coef_offset = 1 + coef_dim1 * (1 + coef_dim2);
    coef -= coef_offset;
    scrtch_dim1 = *k;
    scrtch_dim2 = *k;
    scrtch_offset = 1 + scrtch_dim1 * (1 + scrtch_dim2);
    scrtch -= scrtch_offset;
    bcoef_dim1 = *n;
    bcoef_offset = 1 + bcoef_dim1;
    bcoef -= bcoef_offset;

    /* Function Body */
    lsofar = 0;
    break__[1] = t[*k];
    i__1 = *n;
    for (left = *k; left <= i__1; ++left) {
/*                                find the next nontrivial knot interval. */
	if (t[left + 1] == t[left]) {
	    goto L50;
	}
	++lsofar;
	break__[lsofar + 1] = t[left + 1];
	if (*k > 1) {
	    goto L9;
	}
	i__2 = *m;
	for (mm = 1; mm <= i__2; ++mm) {
/* L5: */
	    coef[mm + (lsofar * coef_dim2 + 1) * coef_dim1] = bcoef[left + mm 
		    * bcoef_dim1];
	}
	goto L50;
/*        store the k b-spline coeff.s relevant to current knot interval */
/*                             in  scrtch(.,1) . */
L9:
	i__2 = *k;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    i__3 = *m;
	    for (mm = 1; mm <= i__3; ++mm) {
/* L10: */
		scrtch[i__ + (mm * scrtch_dim2 + 1) * scrtch_dim1] = bcoef[
			left - *k + i__ + mm * bcoef_dim1];
	    }
	}

/*        for j=1,...,k-1, compute the  k-j  b-spline coeff.s relevant to */
/*        current knot interval for the j-th derivative by differencing */
/*        those for the (j-1)st derivative, and store in scrtch(.,j+1) . */
	i__3 = *k;
	for (jp1 = 2; jp1 <= i__3; ++jp1) {
	    j = jp1 - 1;
	    kmj = *k - j;
	    fkmj = (doublereal) kmj;
	    i__2 = kmj;
	    for (i__ = 1; i__ <= i__2; ++i__) {
		diff = (t[left + i__] - t[left + i__ - kmj]) / fkmj;
		if (diff <= 0.) {
		    goto L20;
		}
		i__4 = *m;
		for (mm = 1; mm <= i__4; ++mm) {
/* L15: */
		    scrtch[i__ + (jp1 + mm * scrtch_dim2) * scrtch_dim1] = (
			    scrtch[i__ + 1 + (j + mm * scrtch_dim2) * 
			    scrtch_dim1] - scrtch[i__ + (j + mm * scrtch_dim2)
			     * scrtch_dim1]) / diff;
		}
L20:
		;
	    }
	}

/*        for  j = 0, ..., k-1, find the values at  t(left)  of the  j+1 */
/*        b-splines of order  j+1  whose support contains the current */
/*        knot interval from those of order  j  (in  biatx ), then comb- */
/*        ine with the b-spline coeff.s (in scrtch(.,k-j) ) found earlier */
/*        to compute the (k-j-1)st derivative at  t(left)  of the given */
/*        spline. */
/*           note. if the repeated calls to  bsplvb  are thought to gene- */
/*        rate too much overhead, then replace the first call by */
/*           biatx(1) = 1. */
/*        and the subsequent call by the statement */
/*           j = jp1 - 1 */
/*        followed by a direct copy of the lines */
/*           deltar(j) = t(left+j) - x */
/*                  ...... */
/*           biatx(j+1) = saved */
/*        from  bsplvb . deltal(kmax)  and  deltar(kmax)  would have to */
/*        appear in a dimension statement, of course. */

	bsplvb_(&t[1], &c__1, &c__1, &t[left], &left, biatx);
	i__2 = *m;
	for (mm = 1; mm <= i__2; ++mm) {
/* L25: */
	    coef[mm + (*k + lsofar * coef_dim2) * coef_dim1] = scrtch[(*k + 
		    mm * scrtch_dim2) * scrtch_dim1 + 1];
	}
	i__2 = *k;
	for (jp1 = 2; jp1 <= i__2; ++jp1) {
	    bsplvb_(&t[1], &jp1, &c__2, &t[left], &left, biatx);
	    kmj = *k + 1 - jp1;
	    i__3 = *m;
	    for (mm = 1; mm <= i__3; ++mm) {
		sum = 0.;
		i__4 = jp1;
		for (i__ = 1; i__ <= i__4; ++i__) {
/* L28: */
		    sum = biatx[i__ - 1] * scrtch[i__ + (kmj + mm * 
			    scrtch_dim2) * scrtch_dim1] + sum;
		}
/* L30: */
		coef[mm + (kmj + lsofar * coef_dim2) * coef_dim1] = sum;
	    }
	}
L50:
	;
    }
    return 0;
} /* bspp2d_ */

doublereal SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::bvalue_(doublereal *t, doublereal *bcoef, integer *n, integer *k, 
	doublereal *x, integer *jderiv)
{
    /* System generated locals */
    integer i__1, i__2;
    doublereal ret_val;

    /* Local variables */
    static integer i__, j;
    static doublereal aj[20];
    static integer jc;
    static doublereal dl[20];
    static integer jj;
    static doublereal dr[20];
    static integer km1, imk, kmj, ilo, nmi;
    static doublereal fkmj;
    static integer mflag, jcmin, jcmax, jdrvp1;

/* *********************************************************************72 */

/* c BVALUE evaluates a derivative of a spline from its B-spline representation. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alls  interv */

/* alculates value at  x  of  jderiv-th derivative of spline from b-repr. */
/*  the spline is taken to be continuous from the right, EXCEPT at the */
/*  rightmost knot, where it is taken to be continuous from the left. */

/* ******  i n p u t ****** */
/*  t, bcoef, n, k......forms the b-representation of the spline  f  to */
/*        be evaluated. specifically, */
/*  t.....knot sequence, of length  n+k, assumed nondecreasing. */
/*  bcoef.....b-coefficient sequence, of length  n . */
/*  n.....length of  bcoef  and dimension of spline(k,t), */
/*        a s s u m e d  positive . */
/*  k.....order of the spline . */

/*  w a r n i n g . . .   the restriction  k .le. kmax (=20)  is imposed */
/*        arbitrarily by the dimension statement for  aj, dl, dr  below, */
/*        but is  n o w h e r e  c h e c k e d  for. */

/*  x.....the point at which to evaluate . */
/*  jderiv.....integer giving the order of the derivative to be evaluated */
/*        a s s u m e d  to be zero or positive. */

/* ******  o u t p u t  ****** */
/*  bvalue.....the value of the (jderiv)-th derivative of  f  at  x . */

/* ******  m e t h o d  ****** */
/*     The nontrivial knot interval  (t(i),t(i+1))  containing  x  is lo- */
/*  cated with the aid of  interv . The  k  b-coeffs of  f  relevant for */
/*  this interval are then obtained from  bcoef (or taken to be zero if */
/*  not explicitly available) and are then differenced  jderiv  times to */
/*  obtain the b-coeffs of  (d**jderiv)f  relevant for that interval. */
/*  Precisely, with  j = jderiv, we have from x.(12) of the text that */

/*     (d**j)f  =  sum ( bcoef(.,j)*b(.,k-j,t) ) */

/*  where */
/*                   / bcoef(.),                     ,  j .eq. 0 */
/*                   / */
/*    bcoef(.,j)  =  / bcoef(.,j-1) - bcoef(.-1,j-1) */
/*                   / ----------------------------- ,  j .gt. 0 */
/*                   /    (t(.+k-j) - t(.))/(k-j) */

/*     Then, we use repeatedly the fact that */

/*    sum ( a(.)*b(.,m,t)(x) )  =  sum ( a(.,x)*b(.,m-1,t)(x) ) */
/*  with */
/*                 (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1) */
/*    a(.,x)  =    --------------------------------------- */
/*                 (x - t(.))      + (t(.+m-1) - x) */

/*  to write  (d**j)f(x)  eventually as a linear combination of b-splines */
/*  of order  1 , and the coefficient for  b(i,1,t)(x)  must then be the */
/*  desired number  (d**j)f(x). (see x.(17)-(19) of text). */

    /* Parameter adjustments */
    --bcoef;
    --t;

    /* Function Body */
    ret_val = 0.;
    if (*jderiv >= *k) {
	goto L99;
    }

/*  *** Find  i   s.t.   1 .le. i .lt. n+k   and   t(i) .lt. t(i+1)   and */
/*      t(i) .le. x .lt. t(i+1) . If no such i can be found,  x  lies */
/*      outside the support of  the spline  f , hence  bvalue = 0. */
/*      (The asymmetry in this choice of  i  makes  f  rightcontinuous, except */
/*      at  t(n+k) where it is leftcontinuous.) */
    i__1 = *n + *k;
    interv_(&t[1], &i__1, x, &i__, &mflag);
    if (mflag != 0) {
	goto L99;
    }
/*  *** if k = 1 (and jderiv = 0), bvalue = bcoef(i). */
    km1 = *k - 1;
    if (km1 > 0) {
	goto L1;
    }
    ret_val = bcoef[i__];
    goto L99;

/*  *** store the k b-spline coefficients relevant for the knot interval */
/*     (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j), */
/*     dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable */
/*     from input to zero. set any t.s not obtainable equal to t(1) or */
/*     to t(n+k) appropriately. */
L1:
    jcmin = 1;
    imk = i__ - *k;
    if (imk >= 0) {
	goto L8;
    }
    jcmin = 1 - imk;
    i__1 = i__;
    for (j = 1; j <= i__1; ++j) {
/* L5: */
	dl[j - 1] = *x - t[i__ + 1 - j];
    }
    i__1 = km1;
    for (j = i__; j <= i__1; ++j) {
	aj[*k - j - 1] = 0.;
/* L6: */
	dl[j - 1] = dl[i__ - 1];
    }
    goto L10;
L8:
    i__1 = km1;
    for (j = 1; j <= i__1; ++j) {
/* L9: */
	dl[j - 1] = *x - t[i__ + 1 - j];
    }

L10:
    jcmax = *k;
    nmi = *n - i__;
    if (nmi >= 0) {
	goto L18;
    }
    jcmax = *k + nmi;
    i__1 = jcmax;
    for (j = 1; j <= i__1; ++j) {
/* L15: */
	dr[j - 1] = t[i__ + j] - *x;
    }
    i__1 = km1;
    for (j = jcmax; j <= i__1; ++j) {
	aj[j] = 0.;
/* L16: */
	dr[j - 1] = dr[jcmax - 1];
    }
    goto L20;
L18:
    i__1 = km1;
    for (j = 1; j <= i__1; ++j) {
/* L19: */
	dr[j - 1] = t[i__ + j] - *x;
    }

L20:
    i__1 = jcmax;
    for (jc = jcmin; jc <= i__1; ++jc) {
/* L21: */
	aj[jc - 1] = bcoef[imk + jc];
    }

/*               *** difference the coefficients  jderiv  times. */
    if (*jderiv == 0) {
	goto L30;
    }
    i__1 = *jderiv;
    for (j = 1; j <= i__1; ++j) {
	kmj = *k - j;
	fkmj = (doublereal) kmj;
	ilo = kmj;
	i__2 = kmj;
	for (jj = 1; jj <= i__2; ++jj) {
	    aj[jj - 1] = (aj[jj] - aj[jj - 1]) / (dl[ilo - 1] + dr[jj - 1]) * 
		    fkmj;
/* L23: */
	    --ilo;
	}
    }

/*  *** compute value at  x  in (t(i),t(i+1)) of jderiv-th derivative, */
/*     given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv). */
L30:
    if (*jderiv == km1) {
	goto L39;
    }
    jdrvp1 = *jderiv + 1;
    i__2 = km1;
    for (j = jdrvp1; j <= i__2; ++j) {
	kmj = *k - j;
	ilo = kmj;
	i__1 = kmj;
	for (jj = 1; jj <= i__1; ++jj) {
	    aj[jj - 1] = (aj[jj] * dl[ilo - 1] + aj[jj - 1] * dr[jj - 1]) / (
		    dl[ilo - 1] + dr[jj - 1]);
/* L33: */
	    --ilo;
	}
    }
L39:
    ret_val = aj[0];

L99:
    return ret_val;
} /* bvalue_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::chol1d_(doublereal *p, doublereal *v, doublereal *qty, 
	integer *npoint, integer *ncol, doublereal *u, doublereal *qu)
{
    /* System generated locals */
    integer v_dim1, v_offset, i__1;

    /* Local variables */
    static integer i__, npm1, npm2;
    static doublereal prev, twop, ratio, six1mp;

/* *********************************************************************72 */

/* c CHOL1D sets up and solves linear systems needed by SMOOTH. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  from  * a practical guide to splines *  by c. de boor */
/* to be called in  s m o o t h */
/* onstructs the upper three diags. in v(i,j), i=2,,npoint-1, j=1,3, of */
/*  the matrix  6*(1-p)*q-transp.*(d**2)*q + p*r, then computes its */
/*  l*l-transp. decomposition and stores it also in v, then applies */
/*  forward and backsubstitution to the right side q-transp.*y in  qty */
/*  to obtain the solution in  u . */

    /* Parameter adjustments */
    --qu;
    --u;
    --qty;
    v_dim1 = *npoint;
    v_offset = 1 + v_dim1;
    v -= v_offset;

    /* Function Body */
    npm1 = *npoint - 1;
/*     construct 6*(1-p)*q-transp.*(d**2)*q  +  p*r */
    six1mp = (1. - *p) * 6.;
    twop = *p * 2.;
    i__1 = npm1;
    for (i__ = 2; i__ <= i__1; ++i__) {
	v[i__ + v_dim1] = six1mp * v[i__ + v_dim1 * 5] + twop * (v[i__ - 1 + (
		v_dim1 << 2)] + v[i__ + (v_dim1 << 2)]);
	v[i__ + (v_dim1 << 1)] = six1mp * v[i__ + v_dim1 * 6] + *p * v[i__ + (
		v_dim1 << 2)];
/* L2: */
	v[i__ + v_dim1 * 3] = six1mp * v[i__ + v_dim1 * 7];
    }
    npm2 = *npoint - 2;
    if (npm2 >= 2) {
	goto L10;
    }
    u[1] = 0.;
    u[2] = qty[2] / v[v_dim1 + 2];
    u[3] = 0.;
    goto L41;
/*  factorization */
L10:
    i__1 = npm2;
    for (i__ = 2; i__ <= i__1; ++i__) {
	ratio = v[i__ + (v_dim1 << 1)] / v[i__ + v_dim1];
	v[i__ + 1 + v_dim1] -= ratio * v[i__ + (v_dim1 << 1)];
	v[i__ + 1 + (v_dim1 << 1)] -= ratio * v[i__ + v_dim1 * 3];
	v[i__ + (v_dim1 << 1)] = ratio;
	ratio = v[i__ + v_dim1 * 3] / v[i__ + v_dim1];
	v[i__ + 2 + v_dim1] -= ratio * v[i__ + v_dim1 * 3];
/* L20: */
	v[i__ + v_dim1 * 3] = ratio;
    }

/*  forward substitution */
    u[1] = 0.;
    v[v_dim1 * 3 + 1] = 0.;
    u[2] = qty[2];
    i__1 = npm2;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L30: */
	u[i__ + 1] = qty[i__ + 1] - v[i__ + (v_dim1 << 1)] * u[i__] - v[i__ - 
		1 + v_dim1 * 3] * u[i__ - 1];
    }
/*  back substitution */
    u[*npoint] = 0.;
    u[npm1] /= v[npm1 + v_dim1];
    i__ = npm2;
L40:
    u[i__] = u[i__] / v[i__ + v_dim1] - u[i__ + 1] * v[i__ + (v_dim1 << 1)] - 
	    u[i__ + 2] * v[i__ + v_dim1 * 3];
    --i__;
    if (i__ > 1) {
	goto L40;
    }
/*  construct q*u */
L41:
    prev = 0.;
    i__1 = *npoint;
    for (i__ = 2; i__ <= i__1; ++i__) {
	qu[i__] = (u[i__] - u[i__ - 1]) / v[i__ - 1 + (v_dim1 << 2)];
	qu[i__ - 1] = qu[i__] - prev;
/* L50: */
	prev = qu[i__];
    }
    qu[*npoint] = -qu[*npoint];
    return 0;
} /* chol1d_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::colloc_(doublereal *aleft, doublereal *aright, integer *
	lbegin, integer *iorder, integer *ntimes, doublereal *addbrk, 
	doublereal *relerr)
{
    /* Format strings */
    static char fmt_630[] = "(\002 approximation from a space of splines of \
order\002,i3,\002 on \002,i3,\002 intervals,\002/\002 of dimension\002,i4\
,\002.  breakpoints -\002/(5e20.10))";
    static char fmt_635[] = "(\002 after\002,i3,\002 of\002,i3,\002 allowed \
iterations,\002)";
    static char fmt_637[] = "(\002 the pp representation of the approximatio\
n is\002)";
    static char fmt_638[] = "(f9.3,e13.6,10e11.3)";
    static char fmt_699[] = "(\002 **********\002/\002 the assigned dimens\
ion\002,i5,\002 for  coef  is too small.\002)";

    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1, d__2, d__3;
    static doublereal equiv_0[2000];

    /* Local variables */
    static doublereal a[200], b[200];
    static integer i__, j, n;
    static doublereal t[200];
    static integer ii;
    static doublereal dx;
    static integer nt;
    static doublereal err, amax;
    static integer iter, lnew, iflag;
    static doublereal asave[200];
#define bloks (equiv_0)
#define templ (equiv_0)
    static doublereal temps[200];
    static integer nncoef;
    static integer nbloks, integs[300]	/* was [3][100] */;

    /* Fortran I/O blocks */
    static cilist io___110 = { 0, 6, 0, fmt_630, 0 };
    static cilist io___111 = { 0, 6, 0, fmt_635, 0 };
    static cilist io___112 = { 0, 6, 0, fmt_637, 0 };
    static cilist io___114 = { 0, 6, 0, fmt_638, 0 };
    static cilist io___118 = { 0, 6, 0, fmt_699, 0 };


/* *********************************************************************72 */

/* c COLLOC solves an ordinary differential equation by collocation. */

/*  from  * a practical guide to splines *  by c. de boor */
/* hapter xv, example. solution of an ode by collocation. */
/* alls colpnt, difequ(ppvalu(interv)), knots, eqblok(putit(difequ*, */
/*     bsplvd(bsplvb)))), slvblk(various subprograms), bsplpp(bsplvb*), */
/*     newnot */

/* ******  i n p u t  ****** */
/*  aleft, aright  endpoints of interval of approximation */
/*  lbegin   initial number of polynomial pieces in the approximation. */
/*           a uniform breakpoint sequence is chosen. */
/*  iorder   order of polynomial pieces in the approximation */
/*  ntimes   number of passes through  n e w n o t  to be made */
/*  addbrk   the number (possibly fractional) of breaks to be added per */
/*           pass through newnot. e.g., if addbrk = .33334, then a break- */
/*           point will be added at every third pass through newnot. */
/*  relerr   a tolerance. newton iteration is stopped if the difference */
/*           between the b-coeffs of two successive iterates is no more */
/*           than  relerr*(absol.largest b-coefficient). */

/* ******  p r i n t e d   o u t p u t  ****** */
/*     consists of the pp-representation of the approximate solution, */
/*     and of the error at selected points. */

/* ******  m e t h o d  ****** */
/*  the m-th order ordinary differential equation with  m  side condit- */
/*  ions, to be specified in subroutine  d i f e q u , is solved approx- */
/*  imately by collocation. */
/*  the approximation  f  to the solution  g  is pp of order  k+m  with */
/*  l  pieces and  m-1 continuous derivatives.  f  is determined by the */
/*  requirement that it satisfy the d.e. at  k  points per interval (to */
/*  be specified in  c o l p n t ) and the  m  side conditions. */
/*     this usually nonlinear system of equations for  f  is solved by */
/*  newton's method. the resulting linear system for the b-coeffs of an */
/*  iterate is constructed appropriately in  e q b l o k  and then solved */
/*  in  s l v b l k , a program designed to solve  a l m o s t  b l o c k */
/*  d i a g o n a l  linear systems efficiently. */
/*     there is an opportunity to attempt improvement of the breakpoint */
/*  sequence (both in number and location) through use of  n e w n o t . */


    approx_1.kpm = *iorder;
    if (*lbegin * approx_1.kpm > 2000) {
	goto L999;
    }
/*  *** set the various parameters concerning the particular dif.equ. */
/*     including a first approx. in case the de is to be solved by */
/*     iteration ( itermx .gt. 0) . */
    difequ_(&c__1, temps, temps);
/*  *** obtain the  k  collocation points for the standard interval. */
    other_1.k = approx_1.kpm - side_1.m;
    colpnt_(&other_1.k, other_1.rho);
/*  *** the following five statements could be replaced by a read in or- */
/*     der to obtain a specific (nonuniform) spacing of the breakpnts. */
    dx = (*aright - *aleft) / (doublereal) (*lbegin);
    temps[0] = *aleft;
    i__1 = *lbegin;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L4: */
	temps[i__ - 1] = temps[i__ - 2] + dx;
    }
    temps[*lbegin] = *aright;
/*  *** generate, in knots, the required knots t(1),...,t(n+kpm). */
    knots_(temps, lbegin, &approx_1.kpm, t, &n);
    nt = 1;
/*  *** generate the almost block diagonal coefficient matrix  bloks  and */
/*     right side  b  from collocation equations and side conditions. */
/*     then solve via  slvblk , obtaining the b-representation of the ap- */
/*     proximation in  t , a ,  n  , kpm  . */
L10:
    eqblok_(t, &n, &approx_1.kpm, temps, a, bloks, &c__2000, integs, &nbloks, 
	    b);
    slvblk_(bloks, integs, &nbloks, b, (integer*)temps, a, &iflag);
    iter = 1;
    if (other_1.itermx <= 1) {
	goto L30;
    }
/*  *** save b-spline coeff. of current approx. in  asave , then get new */
/*     approx. and compare with old. if coeff. are more than  relerr */
/*     apart (relatively) or if no. of iterations is less than  itermx , */
/*     continue iterating. */
L20:
    bsplpp_(t, a, &n, &approx_1.kpm, templ, approx_1.break__, approx_1.coef, &
	    approx_1.l);
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L25: */
	asave[i__ - 1] = a[i__ - 1];
    }
    eqblok_(t, &n, &approx_1.kpm, temps, a, bloks, &c__2000, integs, &nbloks, 
	    b);
    slvblk_(bloks, integs, &nbloks, b, (integer*)temps, a, &iflag);
    err = 0.;
    amax = 0.;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing MAX */
	d__2 = amax, d__3 = (d__1 = a[i__ - 1], abs(d__1));
	amax = max(d__2,d__3);
/* L26: */
/* Computing MAX */
	d__2 = err, d__3 = (d__1 = a[i__ - 1] - asave[i__ - 1], abs(d__1));
	err = max(d__2,d__3);
    }
    if (err <= *relerr * amax) {
	goto L30;
    }
    ++iter;
    if (iter < other_1.itermx) {
	goto L20;
    }
/*  *** iteration (if any) completed. print out approx. based on current */
/*     breakpoint sequence, then try to improve the sequence. */
L30:
    //s_wsfe(&io___110);
    //do_fio(&c__1, (char *)&approx_1.kpm, (ftnlen)sizeof(integer));
    //do_fio(&c__1, (char *)&approx_1.l, (ftnlen)sizeof(integer));
    //do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
    i__1 = approx_1.l;
    //for (i__ = 2; i__ <= i__1; ++i__) {
	//do_fio(&c__1, (char *)&approx_1.break__[i__ - 1], (ftnlen)sizeof(doublereal));
    //}
    //e_wsfe();
    if (other_1.itermx > 0) {
	//s_wsfe(&io___111);
	//do_fio(&c__1, (char *)&iter, (ftnlen)sizeof(integer));
	//do_fio(&c__1, (char *)&other_1.itermx, (ftnlen)sizeof(integer));
	//e_wsfe();
    }
    bsplpp_(t, a, &n, &approx_1.kpm, templ, approx_1.break__, approx_1.coef, &
	    approx_1.l);
    //s_wsfe(&io___112);
    //e_wsfe();
    i__1 = approx_1.l;
    for (i__ = 1; i__ <= i__1; ++i__) {
	ii = (i__ - 1) * approx_1.kpm;
/* L38: */
	//s_wsfe(&io___114);
	//do_fio(&c__1, (char *)&approx_1.break__[i__ - 1], (ftnlen)sizeof(doublereal));
	i__2 = approx_1.kpm;
	/*for (j = 1; j <= i__2; ++j) {
	    do_fio(&c__1, (char *)&approx_1.coef[ii + j - 1], (ftnlen)sizeof(doublereal));
	}*/
	//e_wsfe();
    }
/*  *** the following call is provided here for possible further analysis */
/*     of the approximation specific to the problem being solved. */
/*     it is, of course, easily omitted. */
    difequ_(&c__4, temps, temps);

    if (nt > *ntimes) {
	return 0;
    }
/*  *** from the pp-rep. of the current approx., obtain in  newnot  a new */
/*     (and possibly better) sequence of breakpoints, adding (on the ave- */
/*     rage)  a d d b r k  breakpoints per pass through newnot. */
    lnew = *lbegin + (integer) ((doublereal) nt * *addbrk);
    if (lnew * approx_1.kpm > 2000) {
	goto L999;
    }
    newnot_(approx_1.break__, approx_1.coef, &approx_1.l, &approx_1.kpm, 
	    temps, &lnew, templ);
    knots_(temps, &lnew, &approx_1.kpm, t, &n);
    ++nt;
    goto L10;
L999:
    nncoef = 2000;
    //s_wsfe(&io___118);
    //do_fio(&c__1, (char *)&nncoef, (ftnlen)sizeof(integer));
    //e_wsfe();
    return 0;
} /* colloc_ */

#undef templ
#undef bloks


/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::colpnt_(integer *k, doublereal *rho)
{
    /* Format strings */
    static char fmt_699[] = "(\002 **********\002/\002 equispaced collocatio\
n points are used since k =\002,i2,\002 is greater than 8.\002)";

    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer j;
    static doublereal fkm1o2;

    /* Fortran I/O blocks */
    static cilist io___119 = { 0, 6, 0, fmt_699, 0 };


/* *********************************************************************72 */

/* c COLPNT supplies collocation points. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  the  k collocation points for the standard interval (-1,1) are sup- */
/*  plied here as the zeros of the legendre polynomial of degree k , */
/*  provided  k .le. 8 . otherwise, uniformly spaced points are given. */

    /* Parameter adjustments */
    --rho;

    /* Function Body */
    if (*k > 8) {
	goto L99;
    }
    switch (*k) {
	case 1:  goto L10;
	case 2:  goto L20;
	case 3:  goto L30;
	case 4:  goto L40;
	case 5:  goto L50;
	case 6:  goto L60;
	case 7:  goto L70;
	case 8:  goto L80;
    }
L10:
    rho[1] = 0.;
    return 0;
L20:
    rho[2] = .577350269189626;
    rho[1] = -rho[2];
    return 0;
L30:
    rho[3] = .774596669241483;
    rho[1] = -rho[3];
    rho[2] = 0.;
    return 0;
L40:
    rho[3] = .339981043584856;
    rho[2] = -rho[3];
    rho[4] = .861136311594053;
    rho[1] = -rho[4];
    return 0;
L50:
    rho[4] = .538469310105683;
    rho[2] = -rho[4];
    rho[5] = .906179845938664;
    rho[1] = -rho[5];
    rho[3] = 0.;
    return 0;
L60:
    rho[4] = .238619186083197;
    rho[3] = -rho[4];
    rho[5] = .661209386466265;
    rho[2] = -rho[5];
    rho[6] = .932469514203152;
    rho[1] = -rho[6];
    return 0;
L70:
    rho[5] = .405845151377397;
    rho[3] = -rho[5];
    rho[6] = .741531185599394;
    rho[2] = -rho[6];
    rho[7] = .949107912342759;
    rho[1] = -rho[7];
    rho[4] = 0.;
    return 0;
L80:
    rho[5] = .18343464249565;
    rho[4] = -rho[5];
    rho[6] = .525532409916329;
    rho[3] = -rho[6];
    rho[7] = .796666477413627;
    rho[2] = -rho[7];
    rho[8] = .960289856497536;
    rho[1] = -rho[8];
    return 0;
/*  if k .gt. 8, use equispaced points, but print warning */
L99:
    //s_wsfe(&io___119);
    //do_fio(&c__1, (char *)&(*k), (ftnlen)sizeof(integer));
    //e_wsfe();
    fkm1o2 = (doublereal) (*k - 1) / 2.;
    i__1 = *k;
    for (j = 1; j <= i__1; ++j) {
/* L100: */
	rho[j] = (doublereal) (j - 1) / fkm1o2 - 1.;
    }
    return 0;
} /* colpnt_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::cubspl_(doublereal *tau, doublereal *c__, integer *n, 
	integer *ibcbeg, integer *ibcend)
{
    /* System generated locals */
    integer i__1;
    doublereal d__1;

    /* Local variables */
    static doublereal g;
    static integer i__, j, l, m;
    static doublereal dtau, divdf1, divdf3;

/* *********************************************************************72 */

/* c CUBSPL defines an interpolatory cubic spline. */

/*  from  * a practical guide to splines *  by c. de boor */
/*     ************************  input  *************************** */
/*     n = number of data points. assumed to be .ge. 2. */
/*     (tau(i), c(1,i), i=1,...,n) = abscissae and ordinates of the */
/*        data points. tau is assumed to be strictly increasing. */
/*     ibcbeg, ibcend = boundary condition indicators, and */
/*     c(2,1), c(2,n) = boundary condition information. specifically, */
/*        ibcbeg = 0  means no boundary condition at tau(1) is given. */
/*           in this case, the not-a-knot condition is used, i.e. the */
/*           jump in the third derivative across tau(2) is forced to */
/*           zero, thus the first and the second cubic polynomial pieces */
/*           are made to coincide.) */
/*        ibcbeg = 1  means that the slope at tau(1) is made to equal */
/*           c(2,1), supplied by input. */
/*        ibcbeg = 2  means that the second derivative at tau(1) is */
/*           made to equal c(2,1), supplied by input. */
/*        ibcend = 0, 1, or 2 has analogous meaning concerning the */
/*           boundary condition at tau(n), with the additional infor- */
/*           mation taken from c(2,n). */
/*     ***********************  output  ************************** */
/*     c(j,i), j=1,...,4; i=1,...,l (= n-1) = the polynomial coefficients */
/*        of the cubic interpolating spline with interior knots (or */
/*        joints) tau(2), ..., tau(n-1). precisely, in the interval */
/*        (tau(i), tau(i+1)), the spline f is given by */
/*           f(x) = c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.) */
/*        where h = x - tau(i). the function program *ppvalu* may be */
/*        used to evaluate f or its derivatives from tau,c, l = n-1, */
/*        and k=4. */

/* ****** a tridiagonal linear system for the unknown slopes s(i) of */
/*  f  at tau(i), i=1,...,n, is generated and then solved by gauss elim- */
/*  ination, with s(i) ending up in c(2,i), all i. */
/*     c(3,.) and c(4,.) are used initially for temporary storage. */
    /* Parameter adjustments */
    c__ -= 5;
    --tau;

    /* Function Body */
    l = *n - 1;
/* ompute first differences of tau sequence and store in c(3,.). also, */
/* ompute first divided difference of data and store in c(4,.). */
    i__1 = *n;
    for (m = 2; m <= i__1; ++m) {
	c__[(m << 2) + 3] = tau[m] - tau[m - 1];
/* L10: */
	c__[(m << 2) + 4] = (c__[(m << 2) + 1] - c__[(m - 1 << 2) + 1]) / c__[
		(m << 2) + 3];
    }
/* onstruct first equation from the boundary condition, of the form */
/*             c(4,1)*s(1) + c(3,1)*s(2) = c(2,1) */
    if ((i__1 = *ibcbeg - 1) < 0) {
	goto L11;
    } else if (i__1 == 0) {
	goto L15;
    } else {
	goto L16;
    }
L11:
    if (*n > 2) {
	goto L12;
    }
/*     no condition at left end and n = 2. */
    c__[8] = 1.;
    c__[7] = 1.;
    c__[6] = c__[12] * 2.;
    goto L25;
/*     not-a-knot condition at left end and n .gt. 2. */
L12:
    c__[8] = c__[15];
    c__[7] = c__[11] + c__[15];
/* Computing 2nd power */
    d__1 = c__[11];
    c__[6] = ((c__[11] + c__[7] * 2.) * c__[12] * c__[15] + d__1 * d__1 * c__[
	    16]) / c__[7];
    goto L19;
/*     slope prescribed at left end. */
L15:
    c__[8] = 1.;
    c__[7] = 0.;
    goto L18;
/*     second derivative prescribed at left end. */
L16:
    c__[8] = 2.;
    c__[7] = 1.;
    c__[6] = c__[12] * 3. - c__[11] / 2. * c__[6];
L18:
    if (*n == 2) {
	goto L25;
    }
/*  if there are interior knots, generate the corresp. equations and car- */
/*  ry out the forward pass of gauss elimination, after which the m-th */
/*  equation reads    c(4,m)*s(m) + c(3,m)*s(m+1) = c(2,m). */
L19:
    i__1 = l;
    for (m = 2; m <= i__1; ++m) {
	g = -c__[(m + 1 << 2) + 3] / c__[(m - 1 << 2) + 4];
	c__[(m << 2) + 2] = g * c__[(m - 1 << 2) + 2] + (c__[(m << 2) + 3] * 
		c__[(m + 1 << 2) + 4] + c__[(m + 1 << 2) + 3] * c__[(m << 2) 
		+ 4]) * 3.;
/* L20: */
	c__[(m << 2) + 4] = g * c__[(m - 1 << 2) + 3] + (c__[(m << 2) + 3] + 
		c__[(m + 1 << 2) + 3]) * 2.;
    }
/* onstruct last equation from the second boundary condition, of the form */
/*           (-g*c(4,n-1))*s(n-1) + c(4,n)*s(n) = c(2,n) */
/*     if slope is prescribed at right end, one can go directly to back- */
/*     substitution, since c array happens to be set up just right for it */
/*     at this point. */
    if ((i__1 = *ibcend - 1) < 0) {
	goto L21;
    } else if (i__1 == 0) {
	goto L30;
    } else {
	goto L24;
    }
L21:
    if (*n == 3 && *ibcbeg == 0) {
	goto L22;
    }
/*     not-a-knot and n .ge. 3, and either n.gt.3 or  also not-a-knot at */
/*     left end point. */
    g = c__[(*n - 1 << 2) + 3] + c__[(*n << 2) + 3];
/* Computing 2nd power */
    d__1 = c__[(*n << 2) + 3];
    c__[(*n << 2) + 2] = ((c__[(*n << 2) + 3] + g * 2.) * c__[(*n << 2) + 4] *
	     c__[(*n - 1 << 2) + 3] + d__1 * d__1 * (c__[(*n - 1 << 2) + 1] - 
	    c__[(*n - 2 << 2) + 1]) / c__[(*n - 1 << 2) + 3]) / g;
    g = -g / c__[(*n - 1 << 2) + 4];
    c__[(*n << 2) + 4] = c__[(*n - 1 << 2) + 3];
    goto L29;
/*     either (n=3 and not-a-knot also at left) or (n=2 and not not-a- */
/*     knot at left end point). */
L22:
    c__[(*n << 2) + 2] = c__[(*n << 2) + 4] * 2.;
    c__[(*n << 2) + 4] = 1.;
    goto L28;
/*     second derivative prescribed at right endpoint. */
L24:
    c__[(*n << 2) + 2] = c__[(*n << 2) + 4] * 3. + c__[(*n << 2) + 3] / 2. * 
	    c__[(*n << 2) + 2];
    c__[(*n << 2) + 4] = 2.;
    goto L28;
L25:
    if ((i__1 = *ibcend - 1) < 0) {
	goto L26;
    } else if (i__1 == 0) {
	goto L30;
    } else {
	goto L24;
    }
L26:
    if (*ibcbeg > 0) {
	goto L22;
    }
/*     not-a-knot at right endpoint and at left endpoint and n = 2. */
    c__[(*n << 2) + 2] = c__[(*n << 2) + 4];
    goto L30;
L28:
    g = -1. / c__[(*n - 1 << 2) + 4];
/* omplete forward pass of gauss elimination. */
L29:
    c__[(*n << 2) + 4] = g * c__[(*n - 1 << 2) + 3] + c__[(*n << 2) + 4];
    c__[(*n << 2) + 2] = (g * c__[(*n - 1 << 2) + 2] + c__[(*n << 2) + 2]) / 
	    c__[(*n << 2) + 4];
/* arry out back substitution */
L30:
    j = l;
L40:
    c__[(j << 2) + 2] = (c__[(j << 2) + 2] - c__[(j << 2) + 3] * c__[(j + 1 <<
	     2) + 2]) / c__[(j << 2) + 4];
    --j;
    if (j > 0) {
	goto L40;
    }
/* ****** generate cubic coefficients in each interval, i.e., the deriv.s */
/*  at its left endpoint, from value and slope at its endpoints. */
    i__1 = *n;
    for (i__ = 2; i__ <= i__1; ++i__) {
	dtau = c__[(i__ << 2) + 3];
	divdf1 = (c__[(i__ << 2) + 1] - c__[(i__ - 1 << 2) + 1]) / dtau;
	divdf3 = c__[(i__ - 1 << 2) + 2] + c__[(i__ << 2) + 2] - divdf1 * 2.;
	c__[(i__ - 1 << 2) + 3] = (divdf1 - c__[(i__ - 1 << 2) + 2] - divdf3) 
		* 2. / dtau;
/* L50: */
	c__[(i__ - 1 << 2) + 4] = divdf3 / dtau * (6. / dtau);
    }
    return 0;
} /* cubspl_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::cwidth_(doublereal *w, doublereal *b, integer *nequ, 
	integer *ncols, integer *integs, integer *nbloks, doublereal *d__, 
	doublereal *x, integer *iflag)
{
    /* System generated locals */
    integer w_dim1, w_offset, i__1, i__2, i__3, i__4;
    doublereal d__1, d__2, d__3;

    /* Local variables */
    static integer i__, j, ii;
    static doublereal sum;
    static integer jmax;
    static doublereal temp;
    static integer lasti;
    static doublereal ratio;
    static integer istar;
    static doublereal awi1od;
    static integer ipvtp1, lastcl;
    static doublereal colmax;
    static integer lasteq, nrowad, icount, nexteq, ipvteq;
    static doublereal rowmax;

/* *********************************************************************72 */

/* c CWIDTH solves an almost block diagonal linear system. */

/*  this program is a variation of the theme in the algorithm bandet1 */
/*  by martin and wilkinson (numer.math. 9(1976)279-307). it solves */
/*  the linear system */
/*                           a*x  =  b */
/*  of  nequ  equations in case  a  is almost block diagonal with all */
/*  blocks having  ncols  columns using no more storage than it takes to */
/*  store the interesting part of  a . such systems occur in the determ- */
/*  ination of the b-spline coefficients of a spline approximation. */

/*                           parameters */
/*  w     on input, a two-dimensional array of size (nequ,ncols) contain- */
/*        ing the interesting part of the almost block diagonal coeffici- */
/*        ent matrix  a (see description and example below). the array */
/*        integs  describes the storage scheme. */
/*        on output, w  contains the upper triangular factor  u  of the */
/*        lu factorization of a possibly permuted version of  a . in par- */
/*        ticular, the determinant of  a  could now be found as */
/*            iflag*w(1,1)*w(2,1)* ... * w(nequ,1)  . */
/*  b     on input, the right side of the linear system, of length  nequ. */
/*        the contents of  b  are changed during execution. */
/*  nequ  number of equations in system */
/*  ncols block width, i.e., number of columns in each block. */
/*  integs  integer array, of size (2,nequ), describing the block struct- */
/*        ure of  a . */
/*           integs(1,i) = no. of rows in block i               =  nrow */
/*           integs(2,i) = no. of elimination steps in block i */
/*                       = overhang over next block             =  last */
/*  nbloks  number of blocks */
/*  d     work array, to contain row sizes . if storage is scarce, the */
/*        array  x  could be used in the calling sequence for  d . */
/*  x     on output, contains computed solution (if iflag .ne. 0), of */
/*        length  nequ . */
/*  iflag  on output, integer */
/*         = (-1)**(no.of interchanges during elimination) */
/*                if  a  is invertible */
/*         =  0   if  a  is singular */

/*        ------  block structure of  a  ------ */
/*     the interesting part of  a  is taken to consist of  nbloks  con- */
/*  secutive blocks, with the i-th block made up of  nrowi = integs(1,i) */
/*  consecutive rows and  ncols  consecutive columns of  a , and with */
/*  the first  lasti = integs(2,i) columns to the left of the next block. */
/*  these blocks are stored consecutively in the workarray  w . */
/*     for example, here is an 11th order matrix and its arrangement in */
/*  the workarray  w . (the interesting entries of  a  are indicated by */
/*  their row and column index modulo 10.) */

/*                  ---   a   ---                          ---   w   --- */

/*                     nrow1=3 */
/*          11 12 13 14                                     11 12 13 14 */
/*          21 22 23 24                                     21 22 23 24 */
/*          31 32 33 34      nrow2=2                        31 32 33 34 */
/*   last1=2      43 44 45 46                               43 44 45 46 */
/*                53 54 55 56         nrow3=3               53 54 55 56 */
/*         last2=3         66 67 68 69                      66 67 68 69 */
/*                         76 77 78 79                      76 77 78 79 */
/*                         86 87 88 89   nrow4=1            86 87 88 89 */
/*                  last3=1   97 98 99 90   nrow5=2         97 98 99 90 */
/*                     last4=1   08 09 00 01                08 09 00 01 */
/*                               18 19 10 11                18 19 10 11 */
/*                        last5=4 */

/*  for this interpretation of  a  as an almost block diagonal matrix, */
/*  we have  nbloks = 5 , and the integs array is */

/*                        i=  1   2   3   4   5 */
/*                  k= */
/*  integs(k,i) =      1      3   2   3   1   2 */
/*                     2      2   3   1   1   4 */

/*       --------  method  -------- */
/*     gauss elimination with scaled partial pivoting is used, but mult- */
/*  ipliers are  n o t  s a v e d  in order to save storage. rather, the */
/*  right side is operated on during elimination. */
/*     the two parameters */
/*                  i p v t e q   and  l a s t e q */
/*  are used to keep track of the action.  ipvteq is the index of the */
/*  variable to be eliminated next, from equations  ipvteq+1,...,lasteq, */
/*  using equation  ipvteq (possibly after an interchange) as the pivot */
/*  equation. the entries in the pivot column are  a l w a y s  in column */
/*  1 of  w . this is accomplished by putting the entries in rows */
/*  ipvteq+1,...,lasteq  revised by the elimination of the  ipvteq-th */
/*  variable one to the left in  w . in this way, the columns of the */
/*  equations in a given block (as stored in  w ) will be aligned with */
/*  those of the next block at the moment when these next equations be- */
/*  come involved in the elimination process. */
/*     thus, for the above example, the first elimination steps proceed */
/*  as follows. */

/*  *11 12 13 14    11 12 13 14    11 12 13 14    11 12 13 14 */
/*  *21 22 23 24   *22 23 24       22 23 24       22 23 24 */
/*  *31 32 33 34   *32 33 34      *33 34          33 34 */
/*   43 44 45 46    43 44 45 46   *43 44 45 46   *44 45 46        etc. */
/*   53 54 55 56    53 54 55 56   *53 54 55 56   *54 55 56 */
/*   66 67 68 69    66 67 68 69    66 67 68 69    66 67 68 69 */
/*        .              .              .              . */

/*     in all other respects, the procedure is standard, including the */
/*  scaled partial pivoting. */

    /* Parameter adjustments */
    --x;
    --d__;
    --b;
    w_dim1 = *nequ;
    w_offset = 1 + w_dim1;
    w -= w_offset;
    integs -= 3;

    /* Function Body */
    *iflag = 1;
    ipvteq = 0;
    lasteq = 0;
/*                                 the i-loop runs over the blocks */
    i__1 = *nbloks;
    for (i__ = 1; i__ <= i__1; ++i__) {

/*        the equations for the current block are added to those current- */
/*        ly involved in the elimination process, by increasing  lasteq */
/*        by  integs(1,i) after the rowsize of these equations has been */
/*        recorded in the array  d . */

	nrowad = integs[(i__ << 1) + 1];
	i__2 = nrowad;
	for (icount = 1; icount <= i__2; ++icount) {
	    nexteq = lasteq + icount;
	    rowmax = 0.;
	    i__3 = *ncols;
	    for (j = 1; j <= i__3; ++j) {
/* L5: */
/* Computing MAX */
		d__2 = rowmax, d__3 = (d__1 = w[nexteq + j * w_dim1], abs(
			d__1));
		rowmax = max(d__2,d__3);
	    }
	    if (rowmax == 0.) {
		goto L999;
	    }
/* L10: */
	    d__[nexteq] = rowmax;
	}
	lasteq += nrowad;

/*        there will be  lasti = integs(2,i)  elimination steps before */
/*        the equations in the next block become involved. further, */
/*        l a s t c l  records the number of columns involved in the cur- */
/*        rent elimination step. it starts equal to  ncols  when a block */
/*        first becomes involved and then drops by one after each elim- */
/*        ination step. */

	lastcl = *ncols;
	lasti = integs[(i__ << 1) + 2];
	i__2 = lasti;
	for (icount = 1; icount <= i__2; ++icount) {
	    ++ipvteq;
	    if (ipvteq < lasteq) {
		goto L11;
	    }
	    if ((d__1 = w[ipvteq + w_dim1], abs(d__1)) + d__[ipvteq] > d__[
		    ipvteq]) {
		goto L50;
	    }
	    goto L999;

/*        determine the smallest  i s t a r  in  (ipvteq,lasteq)  for */
/*        which  abs(w(istar,1))/d(istar)  is as large as possible, and */
/*        interchange equations  ipvteq  and  istar  in case  ipvteq */
/*        .lt. istar . */

L11:
	    colmax = (d__1 = w[ipvteq + w_dim1], abs(d__1)) / d__[ipvteq];
	    istar = ipvteq;
	    ipvtp1 = ipvteq + 1;
	    i__3 = lasteq;
	    for (ii = ipvtp1; ii <= i__3; ++ii) {
		awi1od = (d__1 = w[ii + w_dim1], abs(d__1)) / d__[ii];
		if (awi1od <= colmax) {
		    goto L13;
		}
		colmax = awi1od;
		istar = ii;
L13:
		;
	    }
	    if ((d__1 = w[istar + w_dim1], abs(d__1)) + d__[istar] == d__[
		    istar]) {
		goto L999;
	    }
	    if (istar == ipvteq) {
		goto L16;
	    }
	    *iflag = -(*iflag);
	    temp = d__[istar];
	    d__[istar] = d__[ipvteq];
	    d__[ipvteq] = temp;
	    temp = b[istar];
	    b[istar] = b[ipvteq];
	    b[ipvteq] = temp;
	    i__3 = lastcl;
	    for (j = 1; j <= i__3; ++j) {
		temp = w[istar + j * w_dim1];
		w[istar + j * w_dim1] = w[ipvteq + j * w_dim1];
/* L14: */
		w[ipvteq + j * w_dim1] = temp;
	    }

/*        subtract the appropriate multiple of equation  ipvteq  from */
/*        equations  ipvteq+1,...,lasteq to make the coefficient of the */
/*        ipvteq-th unknown (presently in column 1 of  w ) zero, but */
/*        store the new coefficients in  w  one to the left from the old. */

L16:
	    i__3 = lasteq;
	    for (ii = ipvtp1; ii <= i__3; ++ii) {
		ratio = w[ii + w_dim1] / w[ipvteq + w_dim1];
		i__4 = lastcl;
		for (j = 2; j <= i__4; ++j) {
/* L18: */
		    w[ii + (j - 1) * w_dim1] = w[ii + j * w_dim1] - ratio * w[
			    ipvteq + j * w_dim1];
		}
		w[ii + lastcl * w_dim1] = 0.;
/* L20: */
		b[ii] -= ratio * b[ipvteq];
	    }
/* L30: */
	    --lastcl;
	}
L50:
	;
    }

/*  at this point,  w  and  b  contain an upper triangular linear system */
/*  equivalent to the original one, with  w(i,j) containing entry */
/*  (i, i-1+j ) of the coefficient matrix. solve this system by backsub- */
/*  stitution, taking into account its block structure. */

/*                          i-loop over the blocks, in reverse order */
    i__ = *nbloks;
L59:
    lasti = integs[(i__ << 1) + 2];
    jmax = *ncols - lasti;
    i__1 = lasti;
    for (icount = 1; icount <= i__1; ++icount) {
	sum = 0.;
	if (jmax == 0) {
	    goto L61;
	}
	i__2 = jmax;
	for (j = 1; j <= i__2; ++j) {
/* L60: */
	    sum += x[ipvteq + j] * w[ipvteq + (j + 1) * w_dim1];
	}
L61:
	x[ipvteq] = (b[ipvteq] - sum) / w[ipvteq + w_dim1];
	++jmax;
/* L70: */
	--ipvteq;
    }
    --i__;
    if (i__ > 0) {
	goto L59;
    }
    return 0;
L999:
    *iflag = 0;
    return 0;
} /* cwidth_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::difequ_(integer *mode, doublereal *xx, doublereal *v)
{
    /* Format strings */
    static char fmt_499[] = "(\002 carrier,s nonlinear perturb. problem\002)";
    static char fmt_610[] = "(\002 eps \002,e20.10)";
    static char fmt_640[] = "(\002 x, g(x)  and  g(x)-f(x)  at selected poin\
ts\002)";
    static char fmt_641[] = "(3e20.10)";

    /* System generated locals */
    integer i__1;
    doublereal d__1, d__2;

    /* Local variables */
    static integer i__;
    static doublereal x, un, ep1, ep2, eps, error, s2ovep, factor;
    static doublereal solutn;

    /* Fortran I/O blocks */
    static cilist io___149 = { 0, 6, 0, fmt_499, 0 };
    static cilist io___151 = { 0, 6, 0, fmt_610, 0 };
    static cilist io___156 = { 0, 6, 0, fmt_640, 0 };
    static cilist io___162 = { 0, 6, 0, fmt_641, 0 };


/* *********************************************************************72 */

/* c DIFEQU returns information about a differential equation. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alls  ppvalu(interv) */
/*  to be called by   c o l l o c ,  p u t i t */
/*  information about the differential equation is dispensed from here */

/* ******  i n p u t  ****** */
/*  mode  an integer indicating the task to be performed. */
/*        = 1   initialization */
/*        = 2   evaluate  de  at  xx */
/*        = 3   specify the next side condition */
/*        = 4   analyze the approximation */
/*  xx a point at which information is wanted */

/* ******  o u t p u t  ****** */
/*  v  depends on the  mode  . see comments below */


/*  this sample of  difequ  is for the example in chapter xv. it is a */
/*  nonlinear second order two point boundary value problem. */

    /* Parameter adjustments */
    --v;

    /* Function Body */
    switch (*mode) {
	case 1:  goto L10;
	case 2:  goto L20;
	case 3:  goto L30;
	case 4:  goto L40;
    }
/*  initialize everything */
/*  i.e. set the order  m  of the dif.equ., the nondecreasing sequence */
/*  xside(i),i=1,...,m, of points at which side cond.s are given and */
/*  anything else necessary. */
L10:
    side_1.m = 2;
    side_1.xside[0] = 0.;
    side_1.xside[1] = 1.;
/*  *** print out heading */
    //s_wsfe(&io___149);
    //e_wsfe();
    eps = .005;
    //s_wsfe(&io___151);
    //do_fio(&c__1, (char *)&eps, (ftnlen)sizeof(doublereal));
    //e_wsfe();
/*  *** set constants used in formula for solution below. */
/* Computing 2nd power */
    d__1 = sqrt(2.) + sqrt(3.);
    factor = d__1 * d__1;
    s2ovep = sqrt(2. / eps);
/*  *** initial guess for newton iteration. un(x) = x*x - 1. */
    approx_1.l = 1;
    approx_1.break__[0] = 0.;
    i__1 = approx_1.kpm;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L16: */
	approx_1.coef[i__ - 1] = 0.;
    }
    approx_1.coef[0] = -1.;
    approx_1.coef[2] = 2.;
    other_1.itermx = 10;
    return 0;

/*  provide value of left side coeff.s and right side at  xx . */
/*  specifically, at  xx  the dif.equ. reads */
/*        v(m+1)d**m + v(m)d**(m-1) + ... + v(1)d**0  =  v(m+2) */
/*  in terms of the quantities v(i),i=1,...,m+2, to be computed here. */
L20:
    v[3] = eps;
    v[2] = 0.;
    un = ppvalu_(approx_1.break__, approx_1.coef, &approx_1.l, &approx_1.kpm, 
	    xx, &c__0);
    v[1] = un * 2.;
/* Computing 2nd power */
    d__1 = un;
    v[4] = d__1 * d__1 + 1.;
    return 0;

/*  provide the  m  side conditions. these conditions are of the form */
/*        v(m+1)d**m + v(m)d**(m-1) + ... + v(1)d**0  =  v(m+2) */
/*  in terms of the quantities v(i),i=1,...,m+2, to be specified here. */
/*  note that v(m+1) = 0  for customary side conditions. */
L30: 
    v[side_1.m + 1] = 0.;
    switch (side_1.iside) {
	case 1:  goto L31;
	case 2:  goto L32;
	case 3:  goto L39;
    }
L31:
    v[2] = 1.;
    v[1] = 0.;
    v[4] = 0.;
    goto L38;
L32:
    v[2] = 0.;
    v[1] = 1.;
    v[4] = 0.;
L38:
    ++side_1.iside;
L39:
    return 0;

/*  calculate the error near the boundary layer at  1. */
L40:
    //s_wsfe(&io___156);
    //e_wsfe();
    x = (float).75;
    for (i__ = 1; i__ <= 9; ++i__) {
	ep1 = exp(s2ovep * (1. - x)) * factor;
	ep2 = exp(s2ovep * (x + 1.)) * factor;
/* Computing 2nd power */
	d__1 = ep1 + 1.;
/* Computing 2nd power */
	d__2 = ep2 + 1.;
	solutn = 12. / (d__1 * d__1) * ep1 + 12. / (d__2 * d__2) * ep2 - 1.;
	error = solutn - ppvalu_(approx_1.break__, approx_1.coef, &approx_1.l,
		 &approx_1.kpm, &x, &c__0);
	//s_wsfe(&io___162);
	//do_fio(&c__1, (char *)&x, (ftnlen)sizeof(doublereal));
	//do_fio(&c__1, (char *)&solutn, (ftnlen)sizeof(doublereal));
	//do_fio(&c__1, (char *)&error, (ftnlen)sizeof(doublereal));
	//e_wsfe();
/* L41: */
	x += .03125;
    }
    return 0;
} /* difequ_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::dtblok_(doublereal *bloks, integer *integs, integer *
	nbloks, integer *ipivot, integer *iflag, doublereal *detsgn, 
	doublereal *detlog)
{
    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1;

    /* Local variables */
    static integer i__, k, ip, last, nrow, index, indexp;

/* *********************************************************************72 */

/* c DTBLOK gets the determinant of an almost block diagonal matrix. */

/*  computes the determinant of an almost block diagonal matrix whose */
/*  plu factorization has been obtained previously in fcblok. */
/*  *** the logarithm of the determinant is computed instead of the */
/*  determinant itself to avoid the danger of overflow or underflow */
/*  inherent in this calculation. */

/* parameters */
/*    bloks, integs, nbloks, ipivot, iflag  are as on return from fcblok. */
/*            in particular, iflag = (-1)**(number of interchanges dur- */
/*            ing factorization) if successful, otherwise iflag = 0. */
/*    detsgn  on output, contains the sign of the determinant. */
/*    detlog  on output, contains the natural logarithm of the determi- */
/*            nant if determinant is not zero. otherwise contains 0. */


    /* Parameter adjustments */
    --bloks;
    integs -= 4;
    --ipivot;

    /* Function Body */
    *detsgn = (doublereal) (*iflag);
    *detlog = 0.;
    if (*iflag == 0) {
	return 0;
    }
    index = 0;
    indexp = 0;
    i__1 = *nbloks;
    for (i__ = 1; i__ <= i__1; ++i__) {
	nrow = integs[i__ * 3 + 1];
	last = integs[i__ * 3 + 3];
	i__2 = last;
	for (k = 1; k <= i__2; ++k) {
	    ip = index + nrow * (k - 1) + ipivot[indexp + k];
	    *detlog += log((d__1 = bloks[ip], abs(d__1)));
/* L1: */
	    *detsgn *= d_sign(&c_b213, &bloks[ip]);
	}
	index = nrow * integs[i__ * 3 + 2] + index;
/* L2: */
	indexp += nrow;
    }
    return 0;
} /* dtblok_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::eqblok_(doublereal *t, integer *n, integer *kpm, 
	doublereal *work1, doublereal *work2, doublereal *bloks, integer *
	lenblk, integer *integs, integer *nbloks, doublereal *b)
{
    /* Format strings */
    static char fmt_699[] = "(\002 **********\002/\002 the assigned dimens\
ion\002,i5,\002 for  bloks  in  colloc  is too small.\002)";

    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer i__, left, nrow, index;
    static integer indexb, isidel;

    /* Fortran I/O blocks */
    static cilist io___176 = { 0, 6, 0, fmt_699, 0 };


/* *********************************************************************72 */

/* c EQBLOK is to be called in COLLOC. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alls putit(difequ,bsplvd(bsplvb)) */
/*  to be called in  c o l l o c */

/* ******  i n p u t  ****** */
/*  t   the knot sequence, of length n+kpm */
/*  n   the dimension of the approximating spline space, i.e., the order */
/*      of the linear system to be constructed. */
/*  kpm = k+m, the order of the approximating spline */
/*  lenblk   the maximum length of the array  bloks  as allowed by the */
/*           dimension statement in  colloc . */

/* ******  w o r k   a r e a s  ****** */
/*  work1    used in  putit, of size (kpm,kpm) */
/*  work2    used in  putit, of size (kpm,m+1) */

/* ******  o u t p u t  ****** */
/*  bloks    the coefficient matrix of the linear system, stored in al- */
/*           most block diagonal form, of size */
/*              kpm*sum(integs(1,i) , i=1,...,nbloks) */
/*  integs   an integer array, of size (3,nbloks), describing the block */
/*           structure. */
/*           integs(1,i)  =  number of rows in block  i */
/*           integs(2,i)  =  number of columns in block  i */
/*           integs(3,i)  =  number of elimination steps which can be */
/*                       carried out in block  i  before pivoting might */
/*                       bring in an equation from the next block. */
/*  nbloks   number of blocks, equals number of polynomial pieces */
/*  b   the right side of the linear system, stored corresponding to the */
/*      almost block diagonal form, of size sum(integs(1,i) , i=1,..., */
/*      nbloks). */

/* ******  m e t h o d  ****** */
/*  each breakpoint interval gives rise to a block in the linear system. */
/*  this block is determined by the  k  colloc.equations in the interval */
/*  with the side conditions (if any) in the interval interspersed ap- */
/*  propriately, and involves the  kpm  b-splines having the interval in */
/*  their support. correspondingly, such a block has  nrow = k + isidel */
/*  rows, with  isidel = number of side conditions in this and the prev- */
/*  ious intervals, and  ncol = kpm  columns. */
/*     further, because the interior knots have multiplicity  k, we can */
/*  carry out (in slvblk)  k  elimination steps in a block before pivot- */
/*  ing might involve an equation from the next block. in the last block, */
/*  of course, all kpm elimination steps will be carried out (in slvblk). */

/*  see the detailed comments in the solveblok package for further in- */
/*  formation about the almost block diagonal form used here. */

    /* Parameter adjustments */
    --b;
    integs -= 4;
    --bloks;
    --work2;
    --work1;
    --t;

    /* Function Body */
    index = 1;
    indexb = 1;
    i__ = 0;
    side_1.iside = 1;
    i__1 = *n;
    i__2 = other_1.k;
    for (left = *kpm; i__2 < 0 ? left >= i__1 : left <= i__1; left += i__2) {
	++i__;
/*        determine integs(.,i) */
	integs[i__ * 3 + 2] = *kpm;
	if (left < *n) {
	    goto L14;
	}
	integs[i__ * 3 + 3] = *kpm;
	isidel = side_1.m;
	goto L16;
L14:
	integs[i__ * 3 + 3] = other_1.k;
/*        at this point,  iside-1  gives the number of side conditions */
/*        incorporated so far. adding to this the side conditions in the */
/*        current interval gives the number  isidel . */
	isidel = side_1.iside - 1;
L15:
	if (isidel == side_1.m) {
	    goto L16;
	}
	if (side_1.xside[isidel] >= t[left + 1]) {
	    goto L16;
	}
	++isidel;
	goto L15;
L16:
	nrow = other_1.k + isidel;
	integs[i__ * 3 + 1] = nrow;
/*        the detailed equations for this block are generated and put */
/*        together in  p u t i t . */
	if (*lenblk < index + nrow * *kpm - 1) {
	    goto L999;
	}
	putit_(&t[1], kpm, &left, &work1[1], &work2[1], &bloks[index], &nrow, 
		&b[indexb]);
	index += nrow * *kpm;
/* L20: */
	indexb += nrow;
    }
    *nbloks = i__;
    return 0;
L999:
    //s_wsfe(&io___176);
    //do_fio(&c__1, (char *)&(*lenblk), (ftnlen)sizeof(integer));
    //e_wsfe();
    //s_stop("", (ftnlen)0);
    return 0;
} /* eqblok_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::evnnot_(doublereal *break__, doublereal *coef, integer *
	l, integer *k, doublereal *brknew, integer *lnew, doublereal *coefg)
{
    /* System generated locals */
    integer coef_dim1, coef_offset, i__1;

    /* Local variables */
    static integer i__;
    static doublereal step;

/* *********************************************************************72 */

/* c EVNNOT is a version of NEWNOT returning uniform knots. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  this is a fake version of  n e w n o t  , of use in example 3 of */
/*  chapter xiv . */
/*  returns  lnew+1  knots in  brknew  which are equidistributed on (a,b) */
/*  = (break(1),break(l+1)) . */

/* ******  i n p u t ****** */
/*  break, coef, l, k.....contains the pp-representation of a certain */
/*        function  f  of order  k . specifically, */
/*        d**(k-1)f(x) = coef(k,i)  for  break(i).le. x .lt.break(i+1) */
/*  lnew.....number of intervals into which the interval (a,b) is to be */
/*        sectioned by the new breakpoint sequence  brknew . */

/* ******  o u t p u t  ****** */
/*  brknew.....array of length  lnew+1  containing the new breakpoint se- */
/*        quence */
/*  coefg.....the coefficient part of the pp-repr.  break, coefg, l, 2 */
/*        for the monotone p.linear function  g  wrto which  brknew  will */
/*        be equidistributed. */

    /* Parameter adjustments */
    --break__;
    coefg -= 3;
    coef_dim1 = *k;
    coef_offset = 1 + coef_dim1;
    coef -= coef_offset;
    --brknew;

    /* Function Body */
    brknew[1] = break__[1];
    brknew[*lnew + 1] = break__[*l + 1];
    step = (break__[*l + 1] - break__[1]) / (doublereal) (*lnew);
    i__1 = *lnew;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L93: */
	brknew[i__] = break__[1] + (doublereal) (i__ - 1) * step;
    }
    return 0;
} /* evnnot_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::factrb_(doublereal *w, integer *ipivot, doublereal *d__, 
	integer *nrow, integer *ncol, integer *last, integer *iflag)
{
    /* System generated locals */
    integer w_dim1, w_offset, i__1, i__2;
    doublereal d__1, d__2, d__3;

    /* Local variables */
    static integer i__, j, k, kp1;
    static doublereal ratio;
    static integer ipivi, ipivk;
    static doublereal awikdi, colmax, rowmax;

/* *********************************************************************72 */

/* c FACTRB constructs a partial PLU factorization. */

/*  adapted from p.132 of 'element.numer.analysis' by conte-de boor */

/*  constructs a partial plu factorization, corresponding to steps 1,..., */
/*   l a s t   in gauss elimination, for the matrix  w  of order */
/*   ( n r o w ,  n c o l ), using pivoting of scaled rows. */

/*  parameters */
/*    w       contains the (nrow,ncol) matrix to be partially factored */
/*            on input, and the partial factorization on output. */
/*    ipivot  an integer array of length nrow containing a record of the */
/*            pivoting strategy used; row ipivot(i) is used during the */
/*            i-th elimination step, i=1,...,last. */
/*    d       a work array of length nrow used to store row sizes */
/*            temporarily. */
/*    nrow    number of rows of w. */
/*    ncol    number of columns of w. */
/*    last    number of elimination steps to be carried out. */
/*    iflag   on output, equals iflag on input times (-1)**(number of */
/*            row interchanges during the factorization process), in */
/*            case no zero pivot was encountered. */
/*            otherwise, iflag = 0 on output. */

/*  initialize ipivot, d */
    /* Parameter adjustments */
    --d__;
    --ipivot;
    w_dim1 = *nrow;
    w_offset = 1 + w_dim1;
    w -= w_offset;

    /* Function Body */
    i__1 = *nrow;
    for (i__ = 1; i__ <= i__1; ++i__) {
	ipivot[i__] = i__;
	rowmax = 0.;
	i__2 = *ncol;
	for (j = 1; j <= i__2; ++j) {
/* L9: */
/* Computing MAX */
	    d__2 = rowmax, d__3 = (d__1 = w[i__ + j * w_dim1], abs(d__1));
	    rowmax = max(d__2,d__3);
	}
	if (rowmax == 0.) {
	    goto L999;
	}
/* L10: */
	d__[i__] = rowmax;
    }
/* gauss elimination with pivoting of scaled rows, loop over k=1,.,last */
    k = 1;
/*        as pivot row for k-th step, pick among the rows not yet used, */
/*        i.e., from rows ipivot(k),...,ipivot(nrow), the one whose k-th */
/*        entry (compared to the row size) is largest. then, if this row */
/*        does not turn out to be row ipivot(k), redefine ipivot(k) ap- */
/*        propriately and record this interchange by changing the sign */
/*        of  i f l a g . */
L11:
    ipivk = ipivot[k];
    if (k == *nrow) {
	goto L21;
    }
    j = k;
    kp1 = k + 1;
    colmax = (d__1 = w[ipivk + k * w_dim1], abs(d__1)) / d__[ipivk];
/*              find the (relatively) largest pivot */
    i__1 = *nrow;
    for (i__ = kp1; i__ <= i__1; ++i__) {
	ipivi = ipivot[i__];
	awikdi = (d__1 = w[ipivi + k * w_dim1], abs(d__1)) / d__[ipivi];
	if (awikdi <= colmax) {
	    goto L15;
	}
	colmax = awikdi;
	j = i__;
L15:
	;
    }
    if (j == k) {
	goto L16;
    }
    ipivk = ipivot[j];
    ipivot[j] = ipivot[k];
    ipivot[k] = ipivk;
    *iflag = -(*iflag);
L16:
/*        if pivot element is too small in absolute value, declare */
/*        matrix to be noninvertible and quit. */
    if ((d__1 = w[ipivk + k * w_dim1], abs(d__1)) + d__[ipivk] <= d__[ipivk]) 
	    {
	goto L999;
    }
/*        otherwise, subtract the appropriate multiple of the pivot */
/*        row from remaining rows, i.e., the rows ipivot(k+1),..., */
/*        ipivot(nrow), to make k-th entry zero. save the multiplier in */
/*        its place. */
    i__1 = *nrow;
    for (i__ = kp1; i__ <= i__1; ++i__) {
	ipivi = ipivot[i__];
	w[ipivi + k * w_dim1] /= w[ipivk + k * w_dim1];
	ratio = -w[ipivi + k * w_dim1];
	i__2 = *ncol;
	for (j = kp1; j <= i__2; ++j) {
/* L20: */
	    w[ipivi + j * w_dim1] = ratio * w[ipivk + j * w_dim1] + w[ipivi + 
		    j * w_dim1];
	}
    }
    k = kp1;
/*        check for having reached the next block. */
    if (k <= *last) {
	goto L11;
    }
    return 0;
/*     if  last  .eq. nrow , check now that pivot element in last row */
/*     is nonzero. */
L21:
    if ((d__1 = w[ipivk + *nrow * w_dim1], abs(d__1)) + d__[ipivk] > d__[
	    ipivk]) {
	return 0;
    }
/*                   singularity flag set */
L999:
    *iflag = 0;
    return 0;
} /* factrb_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::fcblok_(doublereal *bloks, integer *integs, integer *
	nbloks, integer *ipivot, doublereal *scrtch, integer *iflag)
{
    static integer i__, ncol, last, nrow, index;
    static integer indexb;
    static integer indexn;

/* *********************************************************************72 */

/* c FCBLOK supervises the PLU factorization of an almost block diagonal matrix. */

/* alls subroutines  f a c t r b  and  s h i f t b . */

/*   f c b l o k  supervises the plu factorization with pivoting of */
/*  scaled rows of the almost block diagonal matrix stored in the arrays */
/*   b l o k s  and  i n t e g s . */

/*   factrb = subprogram which carries out steps 1,...,last of gauss */
/*            elimination (with pivoting) for an individual block. */
/*   shiftb = subprogram which shifts the remaining rows to the top of */
/*            the next block */

/* parameters */
/*    bloks   an array that initially contains the almost block diagonal */
/*            matrix  a  to be factored, and on return contains the com- */
/*            puted factorization of  a . */
/*    integs  an integer array describing the block structure of  a . */
/*    nbloks  the number of blocks in  a . */
/*    ipivot  an integer array of dimension  sum (integs(1,n) ; n=1, */
/*            ...,nbloks) which, on return, contains the pivoting stra- */
/*            tegy used. */
/*    scrtch  work area required, of length  max (integs(1,n) ; n=1, */
/*            ...,nbloks). */
/*    iflag   output parameter; */
/*            = 0  in case matrix was found to be singular. */
/*            otherwise, */
/*            = (-1)**(number of row interchanges during factorization) */

    /* Parameter adjustments */
    --bloks;
    integs -= 4;
    --ipivot;
    --scrtch;

    /* Function Body */
    *iflag = 1;
    indexb = 1;
    indexn = 1;
    i__ = 1;
/*                        loop over the blocks.  i  is loop index */
L10:
    index = indexn;
    nrow = integs[i__ * 3 + 1];
    ncol = integs[i__ * 3 + 2];
    last = integs[i__ * 3 + 3];
/*        carry out elimination on the i-th block until next block */
/*        enters, i.e., for columns 1,...,last  of i-th block. */
    factrb_(&bloks[index], &ipivot[indexb], &scrtch[1], &nrow, &ncol, &last, 
	    iflag);
/*         check for having reached a singular block or the last block */
    if (*iflag == 0 || i__ == *nbloks) {
	return 0;
    }
    ++i__;
    indexn = nrow * ncol + index;
/*              put the rest of the i-th block onto the next block */
    shiftb_(&bloks[index], &ipivot[indexb], &nrow, &ncol, &last, &bloks[
	    indexn], &integs[i__ * 3 + 1], &integs[i__ * 3 + 2]);
    indexb += nrow;
    goto L10;
} /* fcblok_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::interv_(doublereal *xt, integer *lxt, doublereal *x, 
	integer *left, integer *mflag)
{
    /* Initialized data */

    static integer ilo = 1;

    static integer ihi, istep, middle;

/* *********************************************************************72 */

/* c INTERV brackets a real value in an ascending vector of values. */

/*  from  * a practical guide to splines *  by C. de Boor */
/* omputes  left = max( i :  xt(i) .lt. xt(lxt) .and.  xt(i) .le. x )  . */

/* ******  i n p u t  ****** */
/*  xt.....a real sequence, of length  lxt , assumed to be nondecreasing */
/*  lxt.....number of terms in the sequence  xt . */
/*  x.....the point whose location with respect to the sequence  xt  is */
/*        to be determined. */

/* ******  o u t p u t  ****** */
/*  left, mflag.....both integers, whose value is */

/*   1     -1      if               x .lt.  xt(1) */
/*   i      0      if   xt(i)  .le. x .lt. xt(i+1) */
/*   i      0      if   xt(i)  .lt. x .eq. xt(i+1) .eq. xt(lxt) */
/*   i      1      if   xt(i)  .lt.        xt(i+1) .eq. xt(lxt) .lt. x */

/*        In particular,  mflag = 0  is the 'usual' case.  mflag .ne. 0 */
/*        indicates that  x  lies outside the CLOSED interval */
/*        xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the */
/*        intervals is due to the decision to make all pp functions cont- */
/*        inuous from the right, but, by returning  mflag = 0  even if */
/*        x = xt(lxt), there is the option of having the computed pp function */
/*        continuous from the left at  xt(lxt) . */

/* ******  m e t h o d  ****** */
/*  The program is designed to be efficient in the common situation that */
/*  it is called repeatedly, with  x  taken from an increasing or decrea- */
/*  sing sequence. This will happen, e.g., when a pp function is to be */
/*  graphed. The first guess for  left  is therefore taken to be the val- */
/*  ue returned at the previous call and stored in the  l o c a l  varia- */
/*  ble  ilo . A first check ascertains that  ilo .lt. lxt (this is nec- */
/*  essary since the present call may have nothing to do with the previ- */
/*  ous call). Then, if  xt(ilo) .le. x .lt. xt(ilo+1), we set  left = */
/*  ilo  and are done after just three comparisons. */
/*     Otherwise, we repeatedly double the difference  istep = ihi - ilo */
/*  while also moving  ilo  and  ihi  in the direction of  x , until */
/*                      xt(ilo) .le. x .lt. xt(ihi) , */
/*  after which we use bisection to get, in addition, ilo+1 = ihi . */
/*  left = ilo  is then returned. */

    /* Parameter adjustments */
    --xt;

    /* Function Body */
    ihi = ilo + 1;
    if (ihi < *lxt) {
	goto L20;
    }
    if (*x >= xt[*lxt]) {
	goto L110;
    }
    if (*lxt <= 1) {
	goto L90;
    }
    ilo = *lxt - 1;
    ihi = *lxt;

L20:
    if (*x >= xt[ihi]) {
	goto L40;
    }
    if (*x >= xt[ilo]) {
	goto L100;
    }

/*              **** now x .lt. xt(ilo) . decrease  ilo  to capture  x . */
    istep = 1;
L31:
    ihi = ilo;
    ilo = ihi - istep;
    if (ilo <= 1) {
	goto L35;
    }
    if (*x >= xt[ilo]) {
	goto L50;
    }
    istep <<= 1;
    goto L31;
L35:
    ilo = 1;
    if (*x < xt[1]) {
	goto L90;
    }
    goto L50;
/*              **** now x .ge. xt(ihi) . increase  ihi  to capture  x . */
L40:
    istep = 1;
L41:
    ilo = ihi;
    ihi = ilo + istep;
    if (ihi >= *lxt) {
	goto L45;
    }
    if (*x < xt[ihi]) {
	goto L50;
    }
    istep <<= 1;
    goto L41;
L45:
    if (*x >= xt[*lxt]) {
	goto L110;
    }
    ihi = *lxt;

/*           **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval. */
L50:
    middle = (ilo + ihi) / 2;
    if (middle == ilo) {
	goto L100;
    }
/*     note. it is assumed that middle = ilo in case ihi = ilo+1 . */
    if (*x < xt[middle]) {
	goto L53;
    }
    ilo = middle;
    goto L50;
L53:
    ihi = middle;
    goto L50;
/* **** set output and return. */
L90:
    *mflag = -1;
    *left = 1;
    return 0;
L100:
    *mflag = 0;
    *left = ilo;
    return 0;
L110:
    *mflag = 1;
    if (*x == xt[*lxt]) {
	*mflag = 0;
    }
    *left = *lxt;
L111:
    if (*left == 1) {
	return 0;
    }
    --(*left);
    if (xt[*left] < xt[*lxt]) {
	return 0;
    }
    goto L111;
} /* interv_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::knots_(doublereal *break__, integer *l, integer *kpm, 
	doublereal *t, integer *n)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer j, k, jj, ll, jjj;

/* *********************************************************************72 */

/* c KNOTS is to be called in COLLOC. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  to be called in  c o l l o c . */
/*  constructs from the given breakpoint sequence  b r e a k  the knot */
/*  sequence  t  so that */
/*  spline(k+m,t) = pp(k+m,break) with  m-1  continuous derivatives . */
/*  this means that */
/*  t(1),...,t(n+kpm)  =  break(1) kpm times, then break(2),..., */
/*        break(l) each  k  times, then, finally, break(l+1) kpm times. */

/* ******  i n p u t  ****** */
/*  break(1),...,break(l+1)  breakpoint sequence */
/*  l  number of intervals or pieces */
/*  kpm   = k + m, order of the pp function or spline */

/* ******  o u t p u t  ****** */
/*  t(1),...,t(n+kpm)  the knot sequence. */
/*  n     = l*k + m  =  dimension of  spline(k+m,t). */

    /* Parameter adjustments */
    --t;
    --break__;

    /* Function Body */
    k = *kpm - side_1.m;
    *n = *l * k + side_1.m;
    jj = *n + *kpm;
    jjj = *l + 1;
    i__1 = *kpm;
    for (ll = 1; ll <= i__1; ++ll) {
	t[jj] = break__[jjj];
/* L11: */
	--jj;
    }
    i__1 = *l;
    for (j = 1; j <= i__1; ++j) {
	--jjj;
	i__2 = k;
	for (ll = 1; ll <= i__2; ++ll) {
	    t[jj] = break__[jjj];
/* L12: */
	    --jj;
	}
    }
    i__2 = *kpm;
    for (ll = 1; ll <= i__2; ++ll) {
/* L13: */
	t[ll] = break__[1];
    }
    return 0;
} /* knots_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::l2appr_(doublereal *t, integer *n, integer *k, 
	doublereal *q, doublereal *diag, doublereal *bcoef)
{
    /* System generated locals */
    integer q_dim1, q_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i__, j, jj, ll, mm;
    static doublereal dw;
    static integer left;
    static doublereal biatx[20];
    static integer leftmk;

/* *********************************************************************72 */

/* c L2APPR constructs a weighted L2 spline approximation to given data. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  to be called in main program  l 2 m a i n . */
/* alls subprograms  bsplvb, bchfac/slv */

/* onstructs the (weighted discrete) l2-approximation by splines of order */
/*  k  with knot sequence  t(1), ..., t(n+k)  to given data points */
/*  ( tau(i), gtau(i) ), i=1,...,ntau. the b-spline coefficients */
/*  b c o e f   of the approximating spline are determined from the */
/*  normal equations using cholesky's method. */

/* ******  i n p u t  ****** */
/*  t(1), ..., t(n+k)  the knot sequence */
/*  n.....the dimension of the space of splines of order k with knots t. */
/*  k.....the order */

/*  w a r n i n g  . . .  the restriction   k .le. kmax (= 20)   is impo- */
/*        sed by the arbitrary dimension statement for  biatx  below, but */
/*        is  n o w h e r e   c h e c k e d   for. */

/* ******  w o r k  a r r a y s  ****** */
/*  q....a work array of size (at least) k*n. its first  k  rows are used */
/*       for the  k  lower diagonals of the gramian matrix  c . */
/*  diag.....a work array of length  n  used in bchfac . */

/* ******  i n p u t  via  c o m m o n  /data/  ****** */
/*  ntau.....number of data points */
/*  (tau(i),gtau(i)), i=1,...,ntau     are the  ntau  data points to be */
/*        fitted . */
/*  weight(i), i=1,...,ntau    are the corresponding weights . */

/* ******  o u t p u t  ****** */
/*  bcoef(1), ..., bcoef(n)  the b-spline coeffs. of the l2-appr. */

/* ******  m e t h o d  ****** */
/*  the b-spline coefficients of the l2-appr. are determined as the sol- */
/*  ution of the normal equations */
/*     sum ( (b(i),b(j))*bcoef(j) : j=1,...,n)  = (b(i),g), */
/*                                               i = 1, ..., n . */
/*  here,  b(i)  denotes the i-th b-spline,  g  denotes the function to */
/*  be approximated, and the  i n n e r   p r o d u c t  of two funct- */
/*  ions  f  and  g  is given by */
/*      (f,g)  :=  sum ( f(tau(i))*g(tau(i))*weight(i) : i=1,...,ntau) . */
/*  the arrays  t a u  and  w e i g h t  are given in common block */
/*   d a t a , as is the array  g t a u  containing the sequence */
/*  g(tau(i)), i=1,...,ntau. */
/*  the relevant function values of the b-splines  b(i), i=1,...,n, are */
/*  supplied by the subprogram  b s p l v b . */
/*     the coeff.matrix  c , with */
/*           c(i,j)  :=  (b(i), b(j)), i,j=1,...,n, */
/*  of the normal equations is symmetric and (2*k-1)-banded, therefore */
/*  can be specified by giving its k bands at or below the diagonal. for */
/*  i=1,...,n,  we store */
/*   (b(i),b(j))  =  c(i,j)  in  q(i-j+1,j), j=i,...,min0(i+k-1,n) */
/*  and the right side */
/*   (b(i), g )  in  bcoef(i) . */
/*  since b-spline values are most efficiently generated by finding sim- */
/*  ultaneously the value of  e v e r y  nonzero b-spline at one point, */
/*  the entries of  c  (i.e., of  q ), are generated by computing, for */
/*  each ll, all the terms involving  tau(ll)  simultaneously and adding */
/*  them to all relevant entries. */


    /* Parameter adjustments */
    --bcoef;
    --diag;
    q_dim1 = *k;
    q_offset = 1 + q_dim1;
    q -= q_offset;
    --t;

    /* Function Body */
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	bcoef[j] = 0.;
	i__2 = *k;
	for (i__ = 1; i__ <= i__2; ++i__) {
/* L7: */
	    q[i__ + j * q_dim1] = 0.;
	}
    }
    left = *k;
    leftmk = 0;
    i__2 = i4data_1.ntau;
    for (ll = 1; ll <= i__2; ++ll) {
/*                   locate  l e f t  s.t. tau(ll) in (t(left),t(left+1)) */
L10:
	if (left == *n) {
	    goto L15;
	}
	if (r8data_1.tau[ll - 1] < t[left + 1]) {
	    goto L15;
	}
	++left;
	++leftmk;
	goto L10;
L15:
	bsplvb_(&t[1], k, &c__1, &r8data_1.tau[ll - 1], &left, biatx);
/*        biatx(mm) contains the value of b(left-k+mm) at tau(ll). */
/*        hence, with  dw := biatx(mm)*weight(ll), the number dw*gtau(ll) */
/*        is a summand in the inner product */
/*           (b(left-k+mm), g)  which goes into  bcoef(left-k+mm) */
/*        and the number biatx(jj)*dw is a summand in the inner product */
/*           (b(left-k+jj), b(left-k+mm)), into  q(jj-mm+1,left-k+mm) */
/*        since  (left-k+jj) - (left-k+mm) + 1  =  jj - mm + 1 . */
	i__1 = *k;
	for (mm = 1; mm <= i__1; ++mm) {
	    dw = biatx[mm - 1] * r8data_1.weight[ll - 1];
	    j = leftmk + mm;
	    bcoef[j] = dw * r8data_1.gtau[ll - 1] + bcoef[j];
	    i__ = 1;
	    i__3 = *k;
	    for (jj = mm; jj <= i__3; ++jj) {
		q[i__ + j * q_dim1] = biatx[jj - 1] * dw + q[i__ + j * q_dim1]
			;
/* L20: */
		++i__;
	    }
	}
    }

/*             construct cholesky factorization for  c  in  q , then use */
/*             it to solve the normal equations */
/*                    c*x  =  bcoef */
/*             for  x , and store  x  in  bcoef . */
    bchfac_(&q[q_offset], k, n, &diag[1]);
    bchslv_(&q[q_offset], k, n, &bcoef[1]);
    return 0;
} /* l2appr_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::l2err_(logical *prfun, doublereal *ftau, doublereal *
	error)
{
    /* Format strings */
    static char fmt_615[] = "(///\002 least square error =\002,e20.6/\002 av\
erage error      =\002,e20.6/\002 maximum error      =\002,e20.6//)";
    static char fmt_620[] = "(///14x,\002approximation and scaled error cu\
rve\002/7x,\002data point\002,7x,\002approximation\002,3x,\002deviation x 10\
**\002,i1/(i4,f16.8,f16.8,f17.6))";

    /* System generated locals */
    integer i__1;
    doublereal d__1;

    /* Local variables */
    static integer ie, ll;
    static doublereal err, errl1, errl2, scale, errmax;

    /* Fortran I/O blocks */
    static cilist io___219 = { 0, 6, 0, fmt_615, 0 };
    static cilist io___222 = { 0, 6, 0, fmt_620, 0 };


/* *********************************************************************72 */

/* c L2ERR computes the errors of an L2 approximation. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  this routine is to be called in the main program  l 2 m a i n . */
/* alls subprogram  ppvalu(interv) */
/*  this subroutine computes various errors of the current l2-approxi- */
/*  mation , whose pp-repr. is contained in common block  approx  , */
/*  to the given data contained in common block  data . it prints out */
/*  the average error  e r r l 1 , the l2-error  e r r l 2,  and the */
/*  maximum error  e r r m a x . */

/* ******  i n p u t  ****** */
/*  prfun  a logical value.  if prfun = TRUE, the routine prints out */
/*          the value of the approximation as well as its error at */
/*          every data point. */

/* ******  o u t p u t  ****** */
/*  ftau(1), ..., ftau(ntau),  with  ftau(i)  the approximation  f at */
/*          tau(i), all i. */
/*  error(1), ..., error(ntau),  with  error(i) = scale*(g - f) */
/*          at tau(i), all i. here,  s c a l e  equals  1. in case */
/*          prfun .ne. 'ON' , or the abs.error is greater than 100 some- */
/*          where. otherwise, s c a l e  is such that the maximum of */
/*          abs(error))  over all  i  lies between  10  and  100. this */
/*          makes the printed output more illustrative. */

    /* Parameter adjustments */
    --error;
    --ftau;

    /* Function Body */
    errl1 = 0.;
    errl2 = 0.;
    errmax = 0.;
    i__1 = i4data_1.ntau;
    for (ll = 1; ll <= i__1; ++ll) {
	ftau[ll] = ppvalu_(approx_2.break__, approx_2.coef, &approx_2.l, &
		approx_2.k, &r8data_2.tau[ll - 1], &c__0);
	error[ll] = r8data_2.gtau[ll - 1] - ftau[ll];
	err = (d__1 = error[ll], abs(d__1));
	if (errmax < err) {
	    errmax = err;
	}
	errl1 += err * r8data_2.weight[ll - 1];
/* Computing 2nd power */
	d__1 = err;
	errl2 += d__1 * d__1 * r8data_2.weight[ll - 1];
    }
    errl1 /= r8data_2.totalw;
    errl2 = sqrt(errl2 / r8data_2.totalw);
    //s_wsfe(&io___219);
    //do_fio(&c__1, (char *)&errl2, (ftnlen)sizeof(doublereal));
    //do_fio(&c__1, (char *)&errl1, (ftnlen)sizeof(doublereal));
    //do_fio(&c__1, (char *)&errmax, (ftnlen)sizeof(doublereal));
    //e_wsfe();

/*  Scale error curve and print. */

    if (*prfun) {
	ie = 0;
	scale = 1.;
	if (errmax >= 10.) {
	    goto L18;
	}
	for (ie = 1; ie <= 9; ++ie) {
	    scale *= 10.;
	    if (errmax * scale >= 10.) {
		goto L18;
	    }
	}
L18:
	i__1 = i4data_1.ntau;
	for (ll = 1; ll <= i__1; ++ll) {
/* L19: */
	    error[ll] *= scale;
	}
	//s_wsfe(&io___222);
	//do_fio(&c__1, (char *)&ie, (ftnlen)sizeof(integer));
	i__1 = i4data_1.ntau;
	/*for (ll = 1; ll <= i__1; ++ll) {
	    do_fio(&c__1, (char *)&ll, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&r8data_2.tau[ll - 1], (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&ftau[ll], (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&error[ll], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();*/
    }
    return 0;
} /* l2err_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::l2knts_(doublereal *break__, integer *l, integer *k, 
	doublereal *t, integer *n)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i__, km1;

/* *********************************************************************72 */

/* c L2KNTS converts breakpoints to knots. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  to be called in main program  l 2 m a i n . */
/* onverts the breakpoint sequence  b r e a k   into a corresponding knot */
/*  sequence  t  to allow the repr. of a pp function of order  k  with */
/*  k-2 continuous derivatives as a spline of order  k  with knot */
/*  sequence  t . this means that */
/*  t(1), ..., t(n+k) =  break(1) k times, then break(i), i=2,...,l, each */
/*                       once, then break(l+1) k times . */
/*  therefore,  n = k-1 + l. */

/* ******  i n p u t  ****** */
/*  k     the order */
/*  l     the number of polynomial pieces */
/*  break(1), ...,break(l+1)  the breakpoint sequence */

/* ******  o u t p u t  ****** */
/*  t(1),...,t(n+k)   the knot sequence */
/*  n     the dimension of the corresp. spline space of order  k . */

/*     dimension break(l+1),t(n+k) */
    /* Parameter adjustments */
    --t;
    --break__;

    /* Function Body */
    km1 = *k - 1;
    i__1 = km1;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L5: */
	t[i__] = break__[1];
    }
    i__1 = *l;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L6: */
	t[km1 + i__] = break__[i__];
    }
    *n = km1 + *l;
    i__1 = *k;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L7: */
	t[*n + i__] = break__[*l + 1];
    }
    return 0;
} /* l2knts_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::newnot_(doublereal *break__, doublereal *coef, integer *
	l, integer *k, doublereal *brknew, integer *lnew, doublereal *coefg)
{
    /* Initialized data */

    static integer iprint = 0;

    /* Format strings */
    static char fmt_600[] = "(\002 step =\002,e16.7/(i5,2e16.5))";

    /* System generated locals */
    integer coef_dim1, coef_offset, i__1;
    doublereal d__1;

    /* Local variables */
    static integer i__, j;
    static doublereal dif, step, stepi, difprv, oneovk;

    /* Fortran I/O blocks */
    static cilist io___231 = { 0, 6, 0, fmt_600, 0 };


/* *********************************************************************72 */

/* c NEWNOT returns LNEW+1 knots which are equidistributed on (A,B) */

/*  from  * a practical guide to splines *  by c. de boor */
/*  returns  lnew+1  knots in  brknew  which are equidistributed on (a,b) */
/*  = (break(1),break(l+1)) wrto a certain monotone fctn  g  related to */
/*  the k-th root of the k-th derivative of the pp function  f  whose pp- */
/*  representation is contained in  break, coef, l, k . */

/* ******  i n p u t ****** */
/*  break, coef, l, k.....contains the pp-representation of a certain */
/*        function  f  of order  k . Specifically, */
/*        d**(k-1)f(x) = coef(k,i)  for  break(i).le. x .lt.break(i+1) */
/*  lnew.....number of intervals into which the interval (a,b) is to be */
/*        sectioned by the new breakpoint sequence  brknew . */

/* ******  o u t p u t  ****** */
/*  brknew.....array of length  lnew+1  containing the new breakpoint se- */
/*        quence */
/*  coefg.....the coefficient part of the pp-repr.  break, coefg, l, 2 */
/*        for the monotone p.linear function  g  wrto which  brknew  will */
/*        be equidistributed. */

/* ******  optional  p r i n t e d  o u t p u t  ****** */
/*  coefg.....the pp coeffs of  g  are printed out if  iprint  is set */
/*        .gt. 0  in data statement below. */

/* ******  m e t h o d  ****** */
/*     The k-th derivative of the given pp function  f  does not exist */
/*  (except perhaps as a linear combination of delta functions). Never- */
/*  theless, we construct a p.constant function  h  with breakpoint se- */
/*  quence  break  which is approximately proportional to abs(d**k(f)). */
/*  Specifically, on  (break(i), break(i+1)), */

/*     abs(jump at break(i) of pc)    abs(jump at break(i+1) of pc) */
/* h = --------------------------  +  ---------------------------- */
/*       break(i+1) - break(i-1)         break(i+2) - break(i) */

/*  with  pc  the p.constant (k-1)st derivative of  f . */
/*      Then, the p.linear function  g  is constructed as */

/*    g(x)  =  integral of  h(y)**(1/k)  for  y  from  a  to  x */

/*  and its pp coeffs. stored in  coefg . */
/*     then  brknew  is determined by */

/*        brknew(i)  =  a + g**(-1)((i-1)*step) , i=1,...,lnew+1 */

/*  where  step = g(b)/lnew  and  (a,b) = (break(1),break(l+1)) . */
/*     In the event that  pc = d**(k-1)(f) is constant in  (a,b)  and */
/*  therefore  h = 0 identically,  brknew  is chosen uniformly spaced. */

/*     dimension break(l+1), brknew(lnew+1) */
/* urrent fortran standard makes it impossible to specify the dimension */
/*  of  break   and  brknew  without the introduction of additional */
/*  otherwise superfluous arguments. */
    /* Parameter adjustments */
    --break__;
    coefg -= 3;
    coef_dim1 = *k;
    coef_offset = 1 + coef_dim1;
    coef -= coef_offset;
    --brknew;

    /* Function Body */

    brknew[1] = break__[1];
    brknew[*lnew + 1] = break__[*l + 1];
/*                               if  g  is constant,  brknew  is uniform. */
    if (*l <= 1) {
	goto L90;
    }
/*                        construct the continuous p.linear function  g . */
    oneovk = 1. / (doublereal) (*k);
    coefg[3] = 0.;
    difprv = (d__1 = coef[*k + (coef_dim1 << 1)] - coef[*k + coef_dim1], abs(
	    d__1)) / (break__[3] - break__[1]);
    i__1 = *l;
    for (i__ = 2; i__ <= i__1; ++i__) {
	dif = (d__1 = coef[*k + i__ * coef_dim1] - coef[*k + (i__ - 1) * 
		coef_dim1], abs(d__1)) / (break__[i__ + 1] - break__[i__ - 1])
		;
	d__1 = dif + difprv;
	coefg[(i__ - 1 << 1) + 2] = pow_dd(&d__1, &oneovk);
	coefg[(i__ << 1) + 1] = coefg[(i__ - 1 << 1) + 1] + coefg[(i__ - 1 << 
		1) + 2] * (break__[i__] - break__[i__ - 1]);
/* L10: */
	difprv = dif;
    }
    d__1 = difprv * 2.;
    coefg[(*l << 1) + 2] = pow_dd(&d__1, &oneovk);
/*                                                     step  =  g(b)/lnew */
    step = (coefg[(*l << 1) + 1] + coefg[(*l << 1) + 2] * (break__[*l + 1] - 
	    break__[*l])) / (doublereal) (*lnew);

    if (iprint > 0) {
	//s_wsfe(&io___231);
	//do_fio(&c__1, (char *)&step, (ftnlen)sizeof(doublereal));
	i__1 = *l;
	/*for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&coefg[(i__ << 1) + 1], (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&coefg[(i__ << 1) + 2], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();*/
    }
/*                              if  g  is constant,  brknew  is uniform . */
    if (step <= 0.) {
	goto L90;
    }

/*      For i=2,...,lnew, construct  brknew(i) = a + g**(-1)(stepi), */
/*      with  stepi = (i-1)*step .  this requires inversion of the p.lin- */
/*      ear function  g .  For this,  j  is found so that */
/*         g(break(j)) .le. stepi .le. g(break(j+1)) */
/*      and then */
/*         brknew(i)  =  break(j) + (stepi-g(break(j)))/dg(break(j)) . */
/*      The midpoint is chosen if  dg(break(j)) = 0 . */
    j = 1;
    i__1 = *lnew;
    for (i__ = 2; i__ <= i__1; ++i__) {
	stepi = (doublereal) (i__ - 1) * step;
L21:
	if (j == *l) {
	    goto L27;
	}
	if (stepi <= coefg[(j + 1 << 1) + 1]) {
	    goto L27;
	}
	++j;
	goto L21;
L27:
	if (coefg[(j << 1) + 2] == 0.) {
	    goto L29;
	}
	brknew[i__] = break__[j] + (stepi - coefg[(j << 1) + 1]) / coefg[(j <<
		 1) + 2];
	goto L30;
L29:
	brknew[i__] = (break__[j] + break__[j + 1]) / (float)2.;
L30:
	;
    }
    return 0;

/*                              if  g  is constant,  brknew  is uniform . */
L90:
    step = (break__[*l + 1] - break__[1]) / (doublereal) (*lnew);
    i__1 = *lnew;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L93: */
	brknew[i__] = break__[1] + (doublereal) (i__ - 1) * step;
    }
    return 0;
} /* newnot_ */

doublereal SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::ppvalu_(doublereal *break__, doublereal *coef, integer *l, integer 
	*k, doublereal *x, integer *jderiv)
{
    /* System generated locals */
    integer coef_dim1, coef_offset, i__1;
    doublereal ret_val;

    /* Local variables */
    static doublereal h__;
    static integer i__, m;
    static doublereal fmmjdr;
    static integer ndummy;

/* *********************************************************************72 */

/* c PPVALU evaluates a piecewise polynomial function or its derivative. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alls  interv */
/* alculates value at  x  of  jderiv-th derivative of pp fct from pp-repr */

/* ******  i n p u t  ****** */
/*  break, coef, l, k.....forms the pp-representation of the function  f */
/*        to be evaluated. specifically, the j-th derivative of  f  is */
/*        given by */

/*     (d**j)f(x) = coef(j+1,i) + h*(coef(j+2,i) + h*( ... (coef(k-1,i) + */
/*                             + h*coef(k,i)/(k-j-1))/(k-j-2) ... )/2)/1 */

/*        with  h = x - break(i),  and */

/*       i  =  max( 1 , max( j ,  break(j) .le. x , 1 .le. j .le. l ) ). */

/*  x.....the point at which to evaluate. */
/*  jderiv.....integer giving the order of the derivative to be evaluat- */
/*        ed.  a s s u m e d  to be zero or positive. */

/* ******  o u t p u t  ****** */
/*  ppvalu.....the value of the (jderiv)-th derivative of  f  at  x. */

/* ******  m e t h o d  ****** */
/*     the interval index  i , appropriate for  x , is found through a */
/*  call to  interv . the formula above for the  jderiv-th derivative */
/*  of  f  is then evaluated (by nested multiplication). */

    /* Parameter adjustments */
    --break__;
    coef_dim1 = *k;
    coef_offset = 1 + coef_dim1;
    coef -= coef_offset;

    /* Function Body */
    ret_val = 0.;
    fmmjdr = (doublereal) (*k - *jderiv);
/*              derivatives of order  k  or higher are identically zero. */
    if (fmmjdr <= 0.) {
	goto L99;
    }

/*              find index  i  of largest breakpoint to the left of  x . */
    i__1 = *l + 1;
    interv_(&break__[1], &i__1, x, &i__, &ndummy);

/*      Evaluate  jderiv-th derivative of  i-th polynomial piece at  x . */
    h__ = *x - break__[i__];
    m = *k;
L9:
    ret_val = ret_val / fmmjdr * h__ + coef[m + i__ * coef_dim1];
    --m;
    fmmjdr += -1.;
    if (fmmjdr > 0.) {
	goto L9;
    }
L99:
    return ret_val;
} /* ppvalu_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::putit_(doublereal *t, integer *kpm, integer *left, 
	doublereal *scrtch, doublereal *dbiatx, doublereal *q, integer *nrow, 
	doublereal *b)
{
    /* System generated locals */
    integer dbiatx_dim1, dbiatx_offset, q_dim1, q_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i__, j;
    static doublereal v[20];
    static integer ll;
    static doublereal dx, xm, xx;
    static integer mp1;
    static doublereal sum;
    static integer mode, irow;
    static integer lowrow;

/* *********************************************************************72 */

/* c PUTIT puts together one block of the collocation equation system. */

/*  from  * a practical guide to splines *  by c. de Boor (7 may 92) */
/* alls  bsplvd(bsplvb),difequ(*) */
/*  to be called by  e q b l o k . */

/*  puts together one block of the collocation equation system */

/* ******  i n p u t  ****** */
/*  t     knot sequence, of size  left+kpm (at least) */
/*  kpm   order of spline */
/*  left  integer indicating interval of interest, viz the interval */
/*           (t(left), t(left+1)) */
/*  nrow  number of rows in block to be put together */

/* ******  w o r k  a r e a  ****** */
/*  scrtch   used in bsplvd, of size (kpm,kpm) */
/*  dbiatx   used to contain derivatives of b-splines, of size (kpm,m+1) */
/*           with dbiatx(j,i+1) containing the i-th derivative of the */
/*           j-th b-spline of interest */

/* ******  o u t p u t  ****** */
/*  q  the block, of size (nrow,kpm) */
/*  b  the corresponding piece of the right side, of size (nrow) */

/* ******  m e t h o d  ****** */
/*  the  k  collocation equations for the interval  (t(left),t(left+1)) */
/*  are constructed with the aid of the subroutine  d i f e q u ( 2, ., */
/*  . ) and interspersed (in order) with the side conditions (if any) in */
/*  this interval, using  d i f e q u ( 3, ., . )  for the information. */
/*     the block  q  has  kpm  columns, corresponding to the  kpm  b- */
/*  splines of order  kpm  which have the interval (t(left),t(left+1)) */
/*  in their support. the block's diagonal is part of the diagonal of the */
/*  total system. the first equation in this block not overlapped by the */
/*  preceding block is therefore equation  l o w r o w , with lowrow = 1+ */
/*  number of side conditions in preceding intervals (or blocks). */

    /* Parameter adjustments */
    --t;
    dbiatx_dim1 = *kpm;
    dbiatx_offset = 1 + dbiatx_dim1;
    dbiatx -= dbiatx_offset;
    --scrtch;
    q_dim1 = *nrow;
    q_offset = 1 + q_dim1;
    q -= q_offset;
    --b;

    /* Function Body */
    mp1 = side_1.m + 1;
    i__1 = *kpm;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *nrow;
	for (i__ = 1; i__ <= i__2; ++i__) {
/* L10: */
	    q[i__ + j * q_dim1] = 0.;
	}
    }
    xm = (t[*left + 1] + t[*left]) / 2.;
    dx = (t[*left + 1] - t[*left]) / 2.;

    ll = 1;
    lowrow = side_1.iside;
    i__2 = *nrow;
    for (irow = lowrow; irow <= i__2; ++irow) {
	if (ll > other_1.k) {
	    goto L22;
	}
	mode = 2;
/*        next collocation point is ... */
	xx = xm + dx * other_1.rho[ll - 1];
	++ll;
/*        the corresp.collocation equation is next unless the next side */
/*        condition occurs at a point at, or to the left of, the next */
/*        collocation point. */
	if (side_1.iside > side_1.m) {
	    goto L24;
	}
	if (side_1.xside[side_1.iside - 1] > xx) {
	    goto L24;
	}
	--ll;
L22:
	mode = 3;
	xx = side_1.xside[side_1.iside - 1];
L24:
	difequ_(&mode, &xx, v);
/*        the next equation, a collocation equation (mode = 2) or a side */
/*        condition (mode = 3), reads */
/*   (*)   (v(m+1)*d**m + v(m)*d**(m-1) +...+ v(1)*d**0)f(xx) = v(m+2) */
/*        in terms of the info supplied by  difequ . the corresponding */
/*        equation for the b-coeffs of  f  therefore has the left side of */
/*        (*), evaluated at each of the  kpm  b-splines having  xx  in */
/*        their support, as its  kpm  possibly nonzero coefficients. */
	bsplvd_(&t[1], kpm, &xx, left, &scrtch[1], &dbiatx[dbiatx_offset], &
		mp1);
	i__1 = *kpm;
	for (j = 1; j <= i__1; ++j) {
	    sum = 0.;
	    i__3 = mp1;
	    for (i__ = 1; i__ <= i__3; ++i__) {
/* L25: */
		sum = v[i__ - 1] * dbiatx[j + i__ * dbiatx_dim1] + sum;
	    }
/* L26: */
	    q[irow + j * q_dim1] = sum;
	}
/* L30: */
	b[irow] = v[side_1.m + 1];
    }
    return 0;
} /* putit_ */

doublereal SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::round_(doublereal *x)
{
    /* Initialized data */

    static doublereal flip = -1.;

    /* System generated locals */
    doublereal ret_val;

/* *********************************************************************72 */

/* c ROUND is called to add some noise to data. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alled in example 1  of chapter xiii */

    flip = -flip;
    ret_val = *x + flip * rount_1.size;
    return ret_val;
} /* round_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::sbblok_(doublereal *bloks, integer *integs, integer *
	nbloks, integer *ipivot, doublereal *b, doublereal *x)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i__, j, nbp1, ncol, last, nrow, index;
    static integer indexb, indexx;

/* *********************************************************************72 */

/* c SBBLOK solves a linear system that was factored by FCBLOK. */

/* alls subroutines  s u b f o r  and  s u b b a k . */

/*  supervises the solution (by forward and backward substitution) of */
/*  the linear system  a*x = b  for x, with the plu factorization of  a */
/*  already generated in  f c b l o k .  individual blocks of equations */
/*  are solved via  s u b f o r  and  s u b b a k . */

/* parameters */
/*    bloks, integs, nbloks, ipivot    are as on return from fcblok. */
/*    b       the right side, stored corresponding to the storage of */
/*            the equations. see comments in  s l v b l k  for details. */
/*    x       solution vector */


/*      forward substitution pass */

    /* Parameter adjustments */
    --bloks;
    integs -= 4;
    --ipivot;
    --b;
    --x;

    /* Function Body */
    index = 1;
    indexb = 1;
    indexx = 1;
    i__1 = *nbloks;
    for (i__ = 1; i__ <= i__1; ++i__) {
	nrow = integs[i__ * 3 + 1];
	last = integs[i__ * 3 + 3];
	subfor_(&bloks[index], &ipivot[indexb], &nrow, &last, &b[indexb], &x[
		indexx]);
	index = nrow * integs[i__ * 3 + 2] + index;
	indexb += nrow;
/* L20: */
	indexx += last;
    }

/*     back substitution pass */

    nbp1 = *nbloks + 1;
    i__1 = *nbloks;
    for (j = 1; j <= i__1; ++j) {
	i__ = nbp1 - j;
	nrow = integs[i__ * 3 + 1];
	ncol = integs[i__ * 3 + 2];
	last = integs[i__ * 3 + 3];
	index -= nrow * ncol;
	indexb -= nrow;
	indexx -= last;
/* L30: */
	subbak_(&bloks[index], &ipivot[indexb], &nrow, &ncol, &last, &x[
		indexx]);
    }
    return 0;
} /* sbblok_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::setupq_(doublereal *x, doublereal *dx, doublereal *y, 
	integer *npoint, doublereal *v, doublereal *qty)
{
    /* System generated locals */
    integer v_dim1, v_offset, i__1;
    doublereal d__1, d__2, d__3;

    /* Local variables */
    static integer i__, npm1;
    static doublereal diff, prev;

/* *********************************************************************72 */

/* c SETUPQ is to be called in SMOOTH. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  to be called in  s m o o t h */
/*  put  delx = x(.+1) - x(.)  into  v(.,4), */
/*  put  the three bands of  q-transp*d  into  v(.,1-3), and */
/*  put the three bands of  (d*q)-transp*(d*q)  at and above the diagonal */
/*     into  v(.,5-7) . */
/*     here,  q is  the tridiagonal matrix of order (npoint-2,npoint) */
/*  with general row  1/delx(i) , -1/delx(i) - 1/delx(i+1) , 1/delx(i+1) */
/*  and   d  is the diagonal matrix  with general row  dx(i) . */

    /* Parameter adjustments */
    --qty;
    v_dim1 = *npoint;
    v_offset = 1 + v_dim1;
    v -= v_offset;
    --y;
    --dx;
    --x;

    /* Function Body */
    npm1 = *npoint - 1;
    v[(v_dim1 << 2) + 1] = x[2] - x[1];
    i__1 = npm1;
    for (i__ = 2; i__ <= i__1; ++i__) {
	v[i__ + (v_dim1 << 2)] = x[i__ + 1] - x[i__];
	v[i__ + v_dim1] = dx[i__ - 1] / v[i__ - 1 + (v_dim1 << 2)];
	v[i__ + (v_dim1 << 1)] = -dx[i__] / v[i__ + (v_dim1 << 2)] - dx[i__] /
		 v[i__ - 1 + (v_dim1 << 2)];
/* L11: */
	v[i__ + v_dim1 * 3] = dx[i__ + 1] / v[i__ + (v_dim1 << 2)];
    }
    v[*npoint + v_dim1] = 0.;
    i__1 = npm1;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L12: */
/* Computing 2nd power */
	d__1 = v[i__ + v_dim1];
/* Computing 2nd power */
	d__2 = v[i__ + (v_dim1 << 1)];
/* Computing 2nd power */
	d__3 = v[i__ + v_dim1 * 3];
	v[i__ + v_dim1 * 5] = d__1 * d__1 + d__2 * d__2 + d__3 * d__3;
    }
    if (npm1 < 3) {
	goto L14;
    }
    i__1 = npm1;
    for (i__ = 3; i__ <= i__1; ++i__) {
/* L13: */
	v[i__ - 1 + v_dim1 * 6] = v[i__ - 1 + (v_dim1 << 1)] * v[i__ + v_dim1]
		 + v[i__ - 1 + v_dim1 * 3] * v[i__ + (v_dim1 << 1)];
    }
L14:
    v[npm1 + v_dim1 * 6] = 0.;
    if (npm1 < 4) {
	goto L16;
    }
    i__1 = npm1;
    for (i__ = 4; i__ <= i__1; ++i__) {
/* L15: */
	v[i__ - 2 + v_dim1 * 7] = v[i__ - 2 + v_dim1 * 3] * v[i__ + v_dim1];
    }
L16:
    v[npm1 - 1 + v_dim1 * 7] = 0.;
    v[npm1 + v_dim1 * 7] = 0.;
/* onstruct  q-transp. * y  in  qty. */
    prev = (y[2] - y[1]) / v[(v_dim1 << 2) + 1];
    i__1 = npm1;
    for (i__ = 2; i__ <= i__1; ++i__) {
	diff = (y[i__ + 1] - y[i__]) / v[i__ + (v_dim1 << 2)];
	qty[i__] = diff - prev;
/* L21: */
	prev = diff;
    }
    return 0;
} /* setupq_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::shiftb_(doublereal *ai, integer *ipivot, integer *nrowi, 
	integer *ncoli, integer *last, doublereal *ai1, integer *nrowi1, 
	integer *ncoli1)
{
    /* System generated locals */
    integer ai_dim1, ai_offset, ai1_dim1, ai1_offset, i__1, i__2;

    /* Local variables */
    static integer j, m, ip, jmax, mmax, jmaxp1;

/* *********************************************************************72 */

/* c SHIFTB shifts the rows in the current block. */

/*  shifts the rows in current block, ai, not used as pivot rows, if */
/*  any, i.e., rows ipivot(last+1),...,ipivot(nrowi), onto the first */
/*  mmax = nrow-last rows of the next block, ai1, with column last+j of */
/*  ai  going to column j , j=1,...,jmax=ncoli-last. the remaining col- */
/*  umns of these rows of ai1 are zeroed out. */

/*                             picture */

/*       original situation after         results in a new block i+1 */
/*       last = 2 columns have been       created and ready to be */
/*       done in factrb (assuming no      factored by next factrb call. */
/*       interchanges of rows) */
/*                   1 */
/*              x  x 1x  x  x           x  x  x  x  x */
/*                   1 */
/*              0  x 1x  x  x           0  x  x  x  x */
/*  block i          1                       --------------- */
/*  nrowi = 4   0  0 1x  x  x           0  0 1x  x  x  0  01 */
/*  ncoli = 5        1                       1             1 */
/*  last = 2    0  0 1x  x  x           0  0 1x  x  x  0  01 */
/*  -------------------------------          1             1   new */
/*                   1x  x  x  x  x          1x  x  x  x  x1  block */
/*                   1                       1             1   i+1 */
/*  block i+1        1x  x  x  x  x          1x  x  x  x  x1 */
/*  nrowi1= 5        1                       1             1 */
/*  ncoli1= 5        1x  x  x  x  x          1x  x  x  x  x1 */
/*  -------------------------------          1-------------1 */
/*                   1 */

    /* Parameter adjustments */
    --ipivot;
    ai_dim1 = *nrowi;
    ai_offset = 1 + ai_dim1;
    ai -= ai_offset;
    ai1_dim1 = *nrowi1;
    ai1_offset = 1 + ai1_dim1;
    ai1 -= ai1_offset;

    /* Function Body */
    mmax = *nrowi - *last;
    jmax = *ncoli - *last;
    if (mmax < 1 || jmax < 1) {
	return 0;
    }
/*              put the remainder of block i into ai1 */
    i__1 = mmax;
    for (m = 1; m <= i__1; ++m) {
	ip = ipivot[*last + m];
	i__2 = jmax;
	for (j = 1; j <= i__2; ++j) {
/* L10: */
	    ai1[m + j * ai1_dim1] = ai[ip + (*last + j) * ai_dim1];
	}
    }
    if (jmax == *ncoli1) {
	return 0;
    }
/*              zero out the upper right corner of ai1 */
    jmaxp1 = jmax + 1;
    i__2 = *ncoli1;
    for (j = jmaxp1; j <= i__2; ++j) {
	i__1 = mmax;
	for (m = 1; m <= i__1; ++m) {
/* L20: */
	    ai1[m + j * ai1_dim1] = 0.;
	}
    }
    return 0;
} /* shiftb_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::slvblk_(doublereal *bloks, integer *integs, integer *
	nbloks, doublereal *b, integer *ipivot, doublereal *x, integer *iflag)
{
/* *********************************************************************72 */

/* c SLVBLK solves the almost block diagonal linear system A * x = b. */

/*    this program solves  the  linear system  a*x = b  where a is an */
/*  almost block diagonal matrix.  such almost block diagonal matrices */
/*  arise naturally in piecewise polynomial interpolation or approx- */
/*  imation and in finite element methods for two-point boundary value */
/*  problems.  the plu factorization method is implemented here to take */
/*  advantage of the special structure of such systems for savings in */
/*  computing time and storage requirements. */

/*                  parameters */
/*  bloks   a one-dimenional array, of length */
/*                   sum( integs(1,i)*integs(2,i) ; i = 1,nbloks ) */
/*          on input, contains the blocks of the almost block diagonal */
/*          matrix  a  .  the array integs (see below and the example) */
/*          describes the block structure. */
/*          on output, contains correspondingly the plu factorization */
/*          of  a  (if iflag .ne. 0).  certain of the entries into bloks */
/*          are arbitrary (where the blocks overlap). */
/*  integs  integer array description of the block structure of  a . */
/*            integs(1,i) = no. of rows of block i        =  nrow */
/*            integs(2,i) = no. of colums of block i      =  ncol */
/*            integs(3,i) = no. of elim. steps in block i =  last */
/*                          i  = 1,2,...,nbloks */
/*          the linear system is of order */
/*                n  =  sum ( integs(3,i) , i=1,...,nbloks ), */
/*          but the total number of rows in the blocks is */
/*              nbrows = sum( integs(1,i) ; i = 1,...,nbloks) */
/*  nbloks  number of blocks */
/*  b       right side of the linear system, array of length nbrows. */
/*          certain of the entries are arbitrary, corresponding to */
/*          rows of the blocks which overlap (see block structure and */
/*          the example below). */
/*  ipivot  on output, integer array containing the pivoting sequence */
/*          used. length is nbrows */
/*  x       on output, contains the computed solution (if iflag .ne. 0) */
/*          length is n. */
/*  iflag   on output, integer */
/*            = (-1)**(no. of interchanges during factorization) */
/*                   if  a  is invertible */
/*            = 0    if  a  is singular */

/*                   auxiliary programs */
/*  fcblok (bloks,integs,nbloks,ipivot,scrtch,iflag)  factors the matrix */
/*           a , and is used for this purpose in slvblk. its arguments */
/*          are as in slvblk, except for */
/*              scrtch = a work array of length max(integs(1,i)). */

/*  sbblok (bloks,integs,nbloks,ipivot,b,x)  solves the system a*x = b */
/*          once  a  is factored. this is done automatically by slvblk */
/*          for one right side b, but subsequent solutions may be */
/*          obtained for additional b-vectors. the arguments are all */
/*          as in slvblk. */

/*  dtblok (bloks,integs,nbloks,ipivot,iflag,detsgn,detlog) computes the */
/*          determinant of  a  once slvblk or fcblok has done the fact- */
/*          orization.the first five arguments are as in slvblk. */
/*              detsgn  = sign of the determinant */
/*              detlog  = natural log of the determinant */

/*             ------ block structure of  a  ------ */
/*  the nbloks blocks are stored consecutively in the array  bloks . */
/*  the first block has its (1,1)-entry at bloks(1), and, if the i-th */
/*  block has its (1,1)-entry at bloks(index(i)), then */
/*         index(i+1) = index(i)  +  nrow(i)*ncol(i) . */
/*    the blocks are pieced together to give the interesting part of  a */
/*  as follows.  for i = 1,2,...,nbloks-1, the (1,1)-entry of the next */
/*  block (the (i+1)st block ) corresponds to the (last+1,last+1)-entry */
/*  of the current i-th block.  recall last = integs(3,i) and note that */
/*  this means that */
/*      a. every block starts on the diagonal of  a . */
/*      b. the blocks overlap (usually). the rows of the (i+1)st block */
/*         which are overlapped by the i-th block may be arbitrarily de- */
/*         fined initially. they are overwritten during elimination. */
/*    the right side for the equations in the i-th block are stored cor- */
/*  respondingly as the last entries of a piece of  b  of length  nrow */
/*  (= integs(1,i)) and following immediately in  b  the corresponding */
/*  piece for the right side of the preceding block, with the right side */
/*  for the first block starting at  b(1) . in this, the right side for */
/*  an equation need only be specified once on input, in the first block */
/*  in which the equation appears. */

/*             ------ example and test driver ------ */
/*    the test driver for this package contains an example, a linear */
/*  system of order 11, whose nonzero entries are indicated in the fol- */
/*  lowing schema by their row and column index modulo 10. next to it */
/*  are the contents of the  integs  arrray when the matrix is taken to */
/*  be almost block diagonal with  nbloks = 5, and below it are the five */
/*  blocks. */

/*                      nrow1 = 3, ncol1 = 4 */
/*           11 12 13 14 */
/*           21 22 23 24   nrow2 = 3, ncol2 = 3 */
/*           31 32 33 34 */
/*  last1 = 2      43 44 45 */
/*                 53 54 55            nrow3 = 3, ncol3 = 4 */
/*        last2 = 3         66 67 68 69   nrow4 = 3, ncol4 = 4 */
/*                          76 77 78 79      nrow5 = 4, ncol5 = 4 */
/*                          86 87 88 89 */
/*                 last3 = 1   97 98 99 90 */
/*                    last4 = 1   08 09 00 01 */
/*                                18 19 10 11 */
/*                       last5 = 4 */

/*         actual input to bloks shown by rows of blocks of  a . */
/*      (the ** items are arbitrary, this storage is used by slvblk) */

/*  11 12 13 14  / ** ** **  / 66 67 68 69  / ** ** ** **  / ** ** ** ** */
/*  21 22 23 24 /  43 44 45 /  76 77 78 79 /  ** ** ** ** /  ** ** ** ** */
/*  31 32 33 34/   53 54 55/   86 87 88 89/   97 98 99 90/   08 09 00 01 */
/*                                                           18 19 10 11 */

/*  index = 1      index = 13  index = 22     index = 34     index = 46 */

/*         actual right side values with ** for arbitrary values */
/*  b1 b2 b3 ** b4 b5 b6 b7 b8 ** ** b9 ** ** b10 b11 */

/*  (it would have been more efficient to combine block 3 with block 4) */

/*     in the call to fcblok,  x  is used for temporary storage. */
    /* Parameter adjustments */
    --bloks;
    integs -= 4;
    --b;
    --ipivot;
    --x;

    /* Function Body */
    fcblok_(&bloks[1], &integs[4], nbloks, &ipivot[1], &x[1], iflag);
    if (*iflag == 0) {
	return 0;
    }
    sbblok_(&bloks[1], &integs[4], nbloks, &ipivot[1], &b[1], &x[1]);
    return 0;
} /* slvblk_ */

doublereal SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::smooth_(doublereal *x, doublereal *y, doublereal *dy, integer *
	npoint, doublereal *s, doublereal *v, doublereal *a)
{
    /* System generated locals */
    integer a_dim1, a_offset, v_dim1, v_offset, i__1;
    doublereal ret_val, d__1;

    /* Local variables */
    static integer i__;
    static doublereal p, q, sfq;
    static integer npm1;
    static doublereal oosf, ooss, sixp, utru, prevq;
    static doublereal change, prevsf;

/* *********************************************************************72 */

/* c SMOOTH constructs the cubic smoothing spline to given data. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alls  setupq, chol1d */

/*  constructs the cubic smoothing spline  f  to given data  (x(i),y(i)), */
/*  i=1,...,npoint, which has as small a second derivative as possible */
/*  while */
/*  s(f) = sum( ((y(i)-f(x(i)))/dy(i))**2 , i=1,...,npoint ) .le. s . */

/* ******  i n p u t  ****** */
/*  x(1),...,x(npoint)   data abscissae,  a s s u m e d  to be strictly */
/*        increasing . */
/*  y(1),...,y(npoint)     corresponding data ordinates . */
/*  dy(1),...,dy(npoint)     estimate of uncertainty in data,  a s s u m- */
/*        e d  to be positive . */
/*  npoint.....number of data points,  a s s u m e d  .gt. 1 */
/*  s.....upper bound on the discrete weighted mean square distance of */
/*        the approximation  f  from the data . */

/* ******  w o r k  a r r a y s  ***** */
/*  v.....of size (npoint,7) */
/*  a.....of size (npoint,4) */

/* *****  o u t p u t  ***** */
/*  a(.,1).....contains the sequence of smoothed ordinates . */
/*  a(i,j) = f^(j-1)(x(i)), j=2,3,4, i=1,...,npoint-1 ,  i.e., the */
/*        first three derivatives of the smoothing spline  f  at the */
/*        left end of each of the data intervals . */
/*     w a r n i n g . . .   a  would have to be transposed before it */
/*        could be used in  ppvalu . */

/* ******  m e t h o d  ****** */
/*     The matrices  Q-transp*d  and  Q-transp*D**2*Q  are constructed in */
/*   s e t u p q  from  x  and  dy , as is the vector  qty = Q-transp*y . */
/*  Then, for given  p , the vector  u  is determined in  c h o l 1 d  as */
/*  the solution of the linear system */
/*               (6(1-p)Q-transp*D**2*Q + p*Q)u  = qty  . */
/*  From  u , the smoothing spline  f  (for this choice of smoothing par- */
/*  ameter  p ) is obtained in the sense that */
/*                        f(x(.))  =  y - 6(1-p)D**2*Q*u        and */
/*                      f''(x(.))  =  6*p*u                      . */
/*     The smoothing parameter  p  is found (if possible) so that */
/*                sf(p)  =  s , */
/*  with  sf(p) = s(f) , where  f  is the smoothing spline as it depends */
/*  on  p .  if  s = 0, then p = 1 . if  sf(0) .le. s , then p = 0 . */
/*  Otherwise, the secant method is used to locate an appropriate  p  in */
/*  the open interval  (0,1) . However, straightforward application of */
/*  the secant method, as done in the original version of this program, */
/*  can be very slow and is influenced by the units in which  x  and  y */
/*  are measured, as C. Reinsch has pointed out. Instead, on recommend- */
/*  ation from C. Reinsch, the secant method is applied to the function */
/*           g:q |--> 1/sqrt{sfq(q)} - 1/sqrt{s} , */
/*  with  sfq(q) := sf(q/(1+q)), since  1/sqrt{sfq}  is monotone increasing */
/*  and close to linear for larger  q . One starts at  q = 0  with a */
/*  Newton step, i.e., */
/*                q_0 = 0,  q_1 = -g(0)/g'(0) */
/*  with  g'(0) = -(1/2) sfq(0)^{-3/2} dsfq, where dsfq = -12*u-transp*r*u , */
/*  and  u  as obtained for  p = 0 . Iteration terminates as soon as */
/*   abs(sf - s) .le. .01*s . */


/*     logical test */
/*     parameter (test = .true.) */
/*     integer itercnt */
    /* Parameter adjustments */
    a_dim1 = *npoint;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    v_dim1 = *npoint;
    v_offset = 1 + v_dim1;
    v -= v_offset;
    --dy;
    --y;
    --x;

    /* Function Body */
    setupq_(&x[1], &dy[1], &y[1], npoint, &v[v_offset], &a[(a_dim1 << 2) + 1])
	    ;
    if (*s > 0.) {
	goto L20;
    }
/* L10: */
    p = 1.;
    chol1d_(&p, &v[v_offset], &a[(a_dim1 << 2) + 1], npoint, &c__1, &a[a_dim1 
	    * 3 + 1], &a[a_dim1 + 1]);
    sfq = 0.;
    goto L60;
L20:
    p = 0.;
    chol1d_(&p, &v[v_offset], &a[(a_dim1 << 2) + 1], npoint, &c__1, &a[a_dim1 
	    * 3 + 1], &a[a_dim1 + 1]);
    sfq = 0.;
    i__1 = *npoint;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L21: */
/* Computing 2nd power */
	d__1 = a[i__ + a_dim1] * dy[i__];
	sfq += d__1 * d__1;
    }
    sfq *= 36.;
    if (sfq <= *s) {
	goto L60;
    }
    utru = 0.;
    i__1 = *npoint;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L25: */
/* Computing 2nd power */
	d__1 = a[i__ + a_dim1 * 3];
	utru += v[i__ - 1 + (v_dim1 << 2)] * (a[i__ - 1 + a_dim1 * 3] * (a[
		i__ - 1 + a_dim1 * 3] + a[i__ + a_dim1 * 3]) + d__1 * d__1);
    }
    ooss = 1. / sqrt(*s);
    oosf = 1. / sqrt(sfq);
    q = -(oosf - ooss) * sfq / (utru * 6. * oosf);
/*  secant iteration for the determination of p starts here. */
/*     itercnt = 0 */
    prevq = 0.;
    prevsf = oosf;
L30:
    d__1 = q / (q + 1.);
    chol1d_(&d__1, &v[v_offset], &a[(a_dim1 << 2) + 1], npoint, &c__1, &a[
	    a_dim1 * 3 + 1], &a[a_dim1 + 1]);
    sfq = 0.;
    i__1 = *npoint;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L35: */
/* Computing 2nd power */
	d__1 = a[i__ + a_dim1] * dy[i__];
	sfq += d__1 * d__1;
    }
/* Computing 2nd power */
    d__1 = q + 1.;
    sfq = sfq * 36. / (d__1 * d__1);
    if ((d__1 = sfq - *s, abs(d__1)) <= *s * .01) {
	goto L59;
    }
    oosf = 1. / sqrt(sfq);
    change = (q - prevq) / (oosf - prevsf) * (oosf - ooss);
    prevq = q;
    q -= change;
    prevsf = oosf;
/*     itercnt = itercnt + 1 */
    goto L30;
L59:
    p = q / (q + 1.);
/* orrect value of p has been found. */
/* ompute pol.coefficients from  Q*u (in a(.,1)). */
L60:
    ret_val = sfq;
/*     if (test) then */
/*        print *, 'number of iterations = ', itercnt */
/*     end if */
/*     six1mp = 6./(1.0D+00+q) */
    i__1 = *npoint;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*  61    a(i,1) = y(i) - six1mp*dy(i)**2*a(i,1) */
/* L61: */
/* Computing 2nd power */
	d__1 = dy[i__];
	a[i__ + a_dim1] = y[i__] - (1. - p) * 6. * (d__1 * d__1) * a[i__ + 
		a_dim1];
    }
    sixp = p * 6.;
    i__1 = *npoint;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L62: */
	a[i__ + a_dim1 * 3] *= sixp;
    }
    npm1 = *npoint - 1;
    i__1 = npm1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	a[i__ + (a_dim1 << 2)] = (a[i__ + 1 + a_dim1 * 3] - a[i__ + a_dim1 * 
		3]) / v[i__ + (v_dim1 << 2)];
/* L63: */
	a[i__ + (a_dim1 << 1)] = (a[i__ + 1 + a_dim1] - a[i__ + a_dim1]) / v[
		i__ + (v_dim1 << 2)] - (a[i__ + a_dim1 * 3] + a[i__ + (a_dim1 
		<< 2)] / 3. * v[i__ + (v_dim1 << 2)]) / 2. * v[i__ + (v_dim1 
		<< 2)];
    }
    return ret_val;
} /* smooth_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::spli2d_(doublereal *tau, doublereal *gtau, doublereal *t,
	 integer *n, integer *k, integer *m, doublereal *work, doublereal *q, 
	doublereal *bcoef, integer *iflag)
{
    /* Format strings */
    static char fmt_699[] = "(\002 linear system in  splint  not invertibl\
e\002)";

    /* System generated locals */
    integer bcoef_dim1, bcoef_offset, gtau_dim1, gtau_offset, i__1, i__2;

    /* Local variables */
    static integer i__, j, jj, km1, np1, left, lenq;
    static doublereal taui;
    static integer kpkm2, ilp1mx;

    /* Fortran I/O blocks */
    static cilist io___293 = { 0, 6, 0, fmt_699, 0 };


/* *********************************************************************72 */

/* c SPLI2D produces a interpolatory tensor product spline. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alls bsplvb, banfac/slv */
/*  this is an extended version of  splint , for the use in tensor prod- -------- */
/*  uct interpolation.                                                   -------- */

/*   spli2d  produces the b-spline coeff.s  bcoef(j,.)  of the spline of -------- */
/*   order  k  with knots  t (i), i=1,..., n + k , which takes on the    -------- */
/*   value  gtau (i,j)  at  tau (i), i=1,..., n , j=1,..., m .           -------- */

/* ******  i n p u t  ****** */
/*  tau.....array of length  n , containing data point abscissae. */
/*    a s s u m p t i o n . . .  tau  is strictly increasing */
/*  gtau(.,j)..corresponding array of length  n , containing data point  -------- */
/*        ordinates, j=1,...,m                                           -------- */
/*  t.....knot sequence, of length  n+k */
/*  n.....number of data points and dimension of spline space  s(k,t) */
/*  k.....order of spline */
/*  m.....number of data sets                                            ******** */

/* ******  w o r k   a r e a  ******                                      ******** */
/*  work  a vector of length  n                                          ******** */

/* ******  o u t p u t  ****** */
/*  q.....array of size  (2*k-1)*n , containing the triangular factoriz- */
/*        ation of the coefficient matrix of the linear system for the b- */
/*        coefficients of the spline interpolant. */
/*           the b-coeffs for the interpolant of an additional data set */
/*        (tau(i),htau(i)), i=1,...,n  with the same data abscissae can */
/*        be obtained without going through all the calculations in this */
/*        routine, simply by loading  htau  into  bcoef  and then execut- */
/*        ing the    call banslv ( q, 2*k-1, n, k-1, k-1, bcoef ) */
/*  bcoef.....the b-coefficients of the interpolant, of length  n */
/*  iflag.....an integer indicating success (= 1)  or failure (= 2) */
/*        the linear system to be solved is (theoretically) invertible if */
/*        and only if */
/*              t(i) .lt. tau(i) .lt. tau(i+k),    all i. */
/*        violation of this condition is certain to lead to  iflag = 2 . */

/* ******  m e t h o d  ****** */
/*     the i-th equation of the linear system  a*bcoef = b  for the b-co- */
/*  effs of the interpolant enforces interpolation at  tau(i), i=1,...,n. */
/*  hence,  b(i) = gtau(i), all i, and  a  is a band matrix with  2k-1 */
/*   bands (if it is invertible). */
/*     the matrix  a  is generated row by row and stored, diagonal by di- */
/*  agonal, in the  r o w s  of the array  q , with the main diagonal go- */
/*  ing into row  k .  see comments in the program below. */
/*     the banded system is then solved by a call to  banfac (which con- */
/*  structs the triangular factorization for  a  and stores it again in */
/*   q ), followed by a call to  banslv (which then obtains the solution */
/*   bcoef  by substitution). */
/*     banfac  does no pivoting, since the total positivity of the matrix */
/*  a  makes this unnecessary. */

/*     dimension q(2*k-1,n), t(n+k) */
/* urrent fortran standard makes it impossible to specify precisely the */
/*  dimension of  q  and  t  without the introduction of otherwise super- */
/*  fluous additional arguments. */
    /* Parameter adjustments */
    --t;
    --work;
    --tau;
    bcoef_dim1 = *m;
    bcoef_offset = 1 + bcoef_dim1;
    bcoef -= bcoef_offset;
    gtau_dim1 = *n;
    gtau_offset = 1 + gtau_dim1;
    gtau -= gtau_offset;
    --q;

    /* Function Body */
    np1 = *n + 1;
    km1 = *k - 1;
    kpkm2 = km1 << 1;
    left = *k;
/*                zero out all entries of q */
    lenq = *n * (*k + km1);
    i__1 = lenq;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L5: */
	q[i__] = 0.;
    }

/*  ***   loop over i to construct the  n  interpolation equations */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	taui = tau[i__];
/* Computing MIN */
	i__2 = i__ + *k;
	ilp1mx = min(i__2,np1);
/*        *** find  left  in the closed interval (i,i+k-1) such that */
/*                t(left) .le. tau(i) .lt. t(left+1) */
/*        matrix is singular if this is not possible */
	left = max(left,i__);
	if (taui < t[left]) {
	    goto L998;
	}
L15:
	if (taui < t[left + 1]) {
	    goto L16;
	}
	++left;
	if (left < ilp1mx) {
	    goto L15;
	}
	--left;
	if (taui > t[left + 1]) {
	    goto L998;
	}
/*        *** the i-th equation enforces interpolation at taui, hence */
/*        a(i,j) = b(j,k,t)(taui), all j. only the  k  entries with  j = */
/*        left-k+1,...,left actually might be nonzero. these  k  numbers */
/*        are returned, in  work  (used for temp.storage here), by the   -------- */
/*        following */
L16:
	bsplvb_(&t[1], k, &c__1, &taui, &left, &work[1]);
/*        we therefore want  work(j) = b(left -k+j)(taui) to go into     -------- */
/*        a(i,left-k+j), i.e., into  q(i-(left+j)+2*k,(left+j)-k) since */
/*        a(i+j,j)  is to go into  q(i+k,j), all i,j,  if we consider  q */
/*        as a two-dim. array , with  2*k-1  rows (see comments in */
/*        banfac). in the present program, we treat  q  as an equivalent */
/*        one-dimensional array (because of fortran restrictions on */
/*        dimension statements) . we therefore want  work(j) to go into  -------- */
/*        entry */
/*            i -(left+j) + 2*k + ((left+j) - k-1)*(2*k-1) */
/*                   =  i-left+1 + (left -k)*(2*k-1) + (2*k-2)*j */
/*        of  q . */
	jj = i__ - left + 1 + (left - *k) * (*k + km1);
	i__2 = *k;
	for (j = 1; j <= i__2; ++j) {
	    jj += kpkm2;
/* L30: */
	    q[jj] = work[j];
	}
    }

/*     ***obtain factorization of  a  , stored again in  q. */
    i__2 = *k + km1;
    banfac_(&q[1], &i__2, n, &km1, &km1, iflag);
    switch (*iflag) {
	case 1:  goto L40;
	case 2:  goto L999;
    }
/*     *** solve  a*bcoef = gtau  by backsubstitution */
L40:
    i__2 = *m;
    for (j = 1; j <= i__2; ++j) {
	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
/* L41: */
	    work[i__] = gtau[i__ + j * gtau_dim1];
	}
	i__1 = *k + km1;
	banslv_(&q[1], &i__1, n, &km1, &km1, &work[1]);
	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
/* L50: */
	    bcoef[j + i__ * bcoef_dim1] = work[i__];
	}
    }
    return 0;
L998:
    *iflag = 2;
L999:
    //s_wsfe(&io___293);
    //e_wsfe();
    return 0;
} /* spli2d_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::splint_(doublereal *tau, doublereal *gtau, doublereal *t,
	 integer *n, integer *k, doublereal *q, doublereal *bcoef, integer *
	iflag)
{
    /* Format strings */
    static char fmt_699[] = "(\002 linear system in  splint  not invertibl\
e\002)";

    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer i__, j, jj, km1, np1, left, lenq;
    static doublereal taui;
    static integer kpkm2, ilp1mx;

    /* Fortran I/O blocks */
    static cilist io___304 = { 0, 6, 0, fmt_699, 0 };


/* *********************************************************************72 */

/* c SPLINT produces the B-spline coefficients BCOEF of an interpolating spline. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alls bsplvb, banfac/slv */

/*   splint  produces the b-spline coeff.s  bcoef  of the spline of order */
/*   k  with knots  t(i), i=1,..., n + k , which takes on the value */
/*   gtau(i) at  tau(i), i=1,..., n . */

/* ******  i n p u t  ****** */
/*  tau.....array of length  n , containing data point abscissae. */
/*    a s s u m p t i o n . . .  tau  is strictly increasing */
/*  gtau.....corresponding array of length  n , containing data point or- */
/*        dinates */
/*  t.....knot sequence, of length  n+k */
/*  n.....number of data points and dimension of spline space  s(k,t) */
/*  k.....order of spline */

/* ******  o u t p u t  ****** */
/*  q.....array of size  (2*k-1)*n , containing the triangular factoriz- */
/*        ation of the coefficient matrix of the linear system for the b- */
/*        coefficients of the spline interpolant. */
/*           the b-coeffs for the interpolant of an additional data set */
/*        (tau(i),htau(i)), i=1,...,n  with the same data abscissae can */
/*        be obtained without going through all the calculations in this */
/*        routine, simply by loading  htau  into  bcoef  and then execut- */
/*        ing the    call banslv ( q, 2*k-1, n, k-1, k-1, bcoef ) */
/*  bcoef.....the b-coefficients of the interpolant, of length  n */
/*  iflag.....an integer indicating success (= 1)  or failure (= 2) */
/*        the linear system to be solved is (theoretically) invertible if */
/*        and only if */
/*              t(i) .lt. tau(i) .lt. t(i+k),    all i. */
/*        violation of this condition is certain to lead to  iflag = 2 . */

/* ******  m e t h o d  ****** */
/*     the i-th equation of the linear system  a*bcoef = b  for the b-co- */
/*  effs of the interpolant enforces interpolation at  tau(i), i=1,...,n. */
/*  hence,  b(i) = gtau(i), all i, and  a  is a band matrix with  2k-1 */
/*   bands (if it is invertible). */
/*     the matrix  a  is generated row by row and stored, diagonal by di- */
/*  agonal, in the  r o w s  of the array  q , with the main diagonal go- */
/*  ing into row  k .  see comments in the program below. */
/*     the banded system is then solved by a call to  banfac (which con- */
/*  structs the triangular factorization for  a  and stores it again in */
/*   q ), followed by a call to  banslv (which then obtains the solution */
/*   bcoef  by substitution). */
/*     banfac  does no pivoting, since the total positivity of the matrix */
/*  a  makes this unnecessary. */

/*     dimension q(2*k-1,n), t(n+k) */
/* urrent fortran standard makes it impossible to specify precisely the */
/*  dimension of  q  and  t  without the introduction of otherwise super- */
/*  fluous additional arguments. */
    /* Parameter adjustments */
    --t;
    --bcoef;
    --gtau;
    --tau;
    --q;

    /* Function Body */
    np1 = *n + 1;
    km1 = *k - 1;
    kpkm2 = km1 << 1;
    left = *k;
/*                zero out all entries of q */
    lenq = *n * (*k + km1);
    i__1 = lenq;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L5: */
	q[i__] = 0.;
    }

/*  ***   loop over i to construct the  n  interpolation equations */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	taui = tau[i__];
/* Computing MIN */
	i__2 = i__ + *k;
	ilp1mx = min(i__2,np1);
/*        *** find  left  in the closed interval (i,i+k-1) such that */
/*                t(left) .le. tau(i) .lt. t(left+1) */
/*        matrix is singular if this is not possible */
	left = max(left,i__);
	if (taui < t[left]) {
	    goto L998;
	}
L15:
	if (taui < t[left + 1]) {
	    goto L16;
	}
	++left;
	if (left < ilp1mx) {
	    goto L15;
	}
	--left;
	if (taui > t[left + 1]) {
	    goto L998;
	}
/*        *** the i-th equation enforces interpolation at taui, hence */
/*        a(i,j) = b(j,k,t)(taui), all j. only the  k  entries with  j = */
/*        left-k+1,...,left actually might be nonzero. these  k  numbers */
/*        are returned, in  bcoef (used for temp.storage here), by the */
/*        following */
L16:
	bsplvb_(&t[1], k, &c__1, &taui, &left, &bcoef[1]);
/*        we therefore want  bcoef(j) = b(left-k+j)(taui) to go into */
/*        a(i,left-k+j), i.e., into  q(i-(left+j)+2*k,(left+j)-k) since */
/*        a(i+j,j)  is to go into  q(i+k,j), all i,j,  if we consider  q */
/*        as a two-dim. array , with  2*k-1  rows (see comments in */
/*        banfac). in the present program, we treat  q  as an equivalent */
/*        one-dimensional array (because of fortran restrictions on */
/*        dimension statements) . we therefore want  bcoef(j) to go into */
/*        entry */
/*            i -(left+j) + 2*k + ((left+j) - k-1)*(2*k-1) */
/*                   =  i-left+1 + (left -k)*(2*k-1) + (2*k-2)*j */
/*        of  q . */
	jj = i__ - left + 1 + (left - *k) * (*k + km1);
	i__2 = *k;
	for (j = 1; j <= i__2; ++j) {
	    jj += kpkm2;
/* L30: */
	    q[jj] = bcoef[j];
	}
    }

/*     ***obtain factorization of  a  , stored again in  q. */
    i__2 = *k + km1;
    banfac_(&q[1], &i__2, n, &km1, &km1, iflag);
    switch (*iflag) {
	case 1:  goto L40;
	case 2:  goto L999;
    }
/*     *** solve  a*bcoef = gtau  by backsubstitution */
L40:
    i__2 = *n;
    for (i__ = 1; i__ <= i__2; ++i__) {
/* L41: */
	bcoef[i__] = gtau[i__];
    }
    i__2 = *k + km1;
    banslv_(&q[1], &i__2, n, &km1, &km1, &bcoef[1]);
    return 0;
L998:
    *iflag = 2;
L999:
   //s_wsfe(&io___304);
    //e_wsfe();
    return 0;
} /* splint_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::splopt_(doublereal *tau, integer *n, integer *k, 
	doublereal *scrtch, doublereal *t, integer *iflag)
{
    /* Initialized data */

    static integer newtmx = 10;
    static doublereal tolrte = 1e-6;

    /* Format strings */
    static char fmt_601[] = "(\002 argument n =\002,i4,\002 in  splopt  is l\
ess than k =\002,i3)";
    static char fmt_602[] = "(\002 argument k =\002,i3,\002 in  splopt  is l\
ess than 3\002)";
    static char fmt_644[] = "(\002 c in  splopt  is not invertible\002)";
    static char fmt_653[] = "(\002 no convergence in  splopt  after\002,i3\
,\002 newton steps.\002)";

    /* System generated locals */
    integer i__1, i__2, i__3;
    doublereal d__1, d__2;

    /* Local variables */
    static integer i__, j, l, id, na, nb, nc, nd, ll, nx, km1, kp1;
    static doublereal del;
    static integer kpk, nmk, kpn;
    static doublereal xij, tol, sum;
    static integer left;
    static doublereal sign;
    static integer lenw, kpkm1, nmkm1, index, llmin, llmax;
    static doublereal delmax, floatk;
    static integer leftmk;
    static doublereal signst;
    static integer newton;

    /* Fortran I/O blocks */
    static cilist io___308 = { 0, 6, 0, fmt_601, 0 };
    static cilist io___309 = { 0, 6, 0, fmt_602, 0 };
    static cilist io___339 = { 0, 6, 0, fmt_644, 0 };
    static cilist io___342 = { 0, 6, 0, fmt_653, 0 };


/* *********************************************************************72 */

/* c SPLOPT computes the knots for an optimal recovery scheme. */

/*  from  * a practical guide to splines *  by c. de boor */
/* alls bsplvb, banfac/slv */
/* omputes the knots  t  for the optimal recovery scheme of order  k */
/*  for data at  tau(i), i=1,...,n . */

/* ******  i n p u t  ****** */
/*  tau.....array of length  n , containing the interpolation points. */
/*     a s s u m e d  to be nondecreasing, with tau(i).lt.tau(i+k),all i. */
/*  n.....number of data points . */
/*  k.....order of the optimal recovery scheme to be used . */

/* ******  w o r k  a r r a y  ***** */
/*  scrtch.....array of length  (n-k)(2k+3) + 5k + 3 . the various */
/*        contents are specified in the text below . */

/* ******  o u t p u t  ****** */
/*  iflag.....integer indicating success (=1) or failure (=2) . */
/*     if iflag = 1, then */
/*  t.....array of length  n+k  containing the optimal knots ready for */
/*        use in optimal recovery. specifically,  t(1) = ... = t(k) = */
/*        tau(1)  and  t(n+1) = ... = t(n+k) = tau(n) , while the  n-k */
/*        interior knots  t(k+1), ..., t(n)  are calculated as described */
/*        below under  *method* . */
/*     if iflag = 2, then */
/*        k .lt. 3, or n .lt. k, or a certain linear system was found to */
/*        be singular. */

/* ******  p r i n t e d  o u t p u t  ****** */
/*  a comment will be printed in case  iflag = 2  or newton iterations */
/*  failed to converge in  n e w t m x  iterations . */

/* ******  m e t h o d  ****** */
/*     the (interior) knots  t(k+1), ..., t(n)  are determined by newtons */
/*  method in such a way that the signum function which changes sign at */
/*   t(k+1), ..., t(n)  and nowhere else in  (tau(1),tau(n)) is orthogon- */
/*  al to the spline space  spline( k , tau )  on that interval . */
/*     let  xi(j)  be the current guess for  t(k+j), j=1,...,n-k. then */
/*  the next newton iterate is of the form */
/*              xi(j)  +  (-)**(n-k-j)*x(j)  ,  j=1,...,n-k, */
/*  with  x  the solution of the linear system */
/*                        c*x  =  d  . */
/*  here,  c(i,j) = b(i)(xi(j)), all j, with  b(i)  the i-th b-spline of */
/*  order  k  for the knot sequence  tau , all i, and  d  is the vector */
/*  given by  d(i) = sum( -a(j) , j=i,...,n )*(tau(i+k)-tau(i))/k, all i, */
/*  with  a(i) = sum ( (-)**(n-k-j)*b(i,k+1,tau)(xi(j)) , j=1,...,n-k ) */
/*  for i=1,...,n-1, and  a(n) = -.5 . */
/*     (see chapter  xiii  of text and references there for a derivation) */
/*     the first guess for  t(k+j)  is  (tau(j+1)+...+tau(j+k-1))/(k-1) . */
/*     iteration terminates if  max(abs(x(j))) .lt. t o l  , with */
/*                 t o l  =  t o l r t e *(tau(n)-tau(1))/(n-k) , */
/*  or else after  n e w t m x  iterations , currently, */
/*                 newtmx, tolrte / 10, .000001 */

/*     dimension scrtch((n-k)*(2*k+3)+5*k+3), t(n+k) */
/* urrent fortran standard makes it impossible to specify the precise dim- */
/*  ensions of  scrtch  and  t  without the introduction of otherwise */
/*  superfluous additional arguments . */
    /* Parameter adjustments */
    --tau;
    --scrtch;
    --t;

    /* Function Body */
    nmk = *n - *k;
    if (nmk < 0) {
	goto L1;
    } else if (nmk == 0) {
	goto L56;
    } else {
	goto L2;
    }
L1:
    //s_wsfe(&io___308);
    //do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
    //do_fio(&c__1, (char *)&(*k), (ftnlen)sizeof(integer));
    //e_wsfe();
    goto L999;
L2:
    if (*k > 2) {
	goto L3;
    }
    //s_wsfe(&io___309);
    //do_fio(&c__1, (char *)&(*k), (ftnlen)sizeof(integer));
    //e_wsfe();
    goto L999;
L3:
    nmkm1 = nmk - 1;
    floatk = (doublereal) (*k);
    kpk = *k + *k;
    kp1 = *k + 1;
    km1 = *k - 1;
    kpkm1 = kpk - 1;
    kpn = *k + *n;
    signst = -1.;
    if (nmk > nmk / 2 << 1) {
	signst = 1.;
    }
/*  scrtch(i) = tau-extended(i), i=1,...,n+k+k */
    nx = *n + kpk + 1;
/*  scrtch(i+nx) = xi(i),i=0,...,n-k+1 */
    na = nx + nmk + 1;
/*  scrtch(i+na) = -a(i), i=1,...,n */
    nd = na + *n;
/*  scrtch(i+nd) = x(i) or d(i), i=1,...,n-k */
    nb = nd + nmk;
/*  scrtch(i+nb) = biatx(i),i=1,...,k+1 */
    nc = nb + kp1;
/*  scrtch(i+(j-1)*(2k-1) + nc) = w(i,j) = c(i-k+j,j), i=j-k,...,j+k, */
/*                                                     j=1,...,n-k. */
    lenw = kpkm1 * nmk;
/*  extend  tau  to a knot sequence and store in scrtch. */
    i__1 = *k;
    for (j = 1; j <= i__1; ++j) {
	scrtch[j] = tau[1];
/* L5: */
	scrtch[kpn + j] = tau[*n];
    }
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/* L6: */
	scrtch[*k + j] = tau[j];
    }
/*  first guess for  scrtch (.+nx)  =  xi . */
    scrtch[nx] = tau[1];
    scrtch[nmk + 1 + nx] = tau[*n];
    i__1 = nmk;
    for (j = 1; j <= i__1; ++j) {
	sum = 0.;
	i__2 = km1;
	for (l = 1; l <= i__2; ++l) {
/* L9: */
	    sum += tau[j + l];
	}
/* L10: */
	scrtch[j + nx] = sum / (doublereal) km1;
    }
/*  last entry of  scrtch (.+na)  =  - a  is always ... */
    scrtch[*n + na] = (float).5;
/*  start newton iteration. */
    newton = 1;
    tol = tolrte * (tau[*n] - tau[1]) / (doublereal) nmk;
/*  start newton step */
/* ompute the 2k-1 bands of the matrix c and store in scrtch(.+nc), */
/*  and compute the vector  scrtch(.+na) = -a. */
L20:
    i__1 = lenw;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L21: */
	scrtch[i__ + nc] = 0.;
    }
    i__1 = *n;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L22: */
	scrtch[i__ - 1 + na] = 0.;
    }
    sign = signst;
    left = kp1;
    i__1 = nmk;
    for (j = 1; j <= i__1; ++j) {
	xij = scrtch[j + nx];
L23:
	if (xij < scrtch[left + 1]) {
	    goto L25;
	}
	++left;
	if (left < kpn) {
	    goto L23;
	}
	--left;
L25:
	bsplvb_(&scrtch[1], k, &c__1, &xij, &left, &scrtch[nb + 1]);
/*        the tau sequence in scrtch is preceded by  k  additional knots */
/*        therefore,  scrtch(ll+nb)  now contains  b(left-2k+ll)(xij) */
/*        which is destined for  c(left-2k+ll,j), and therefore for */
/*            w(left-k-j+ll,j)= scrtch(left-k-j+ll + (j-1)*kpkm1 + nc) */
/*        since we store the 2k-1 bands of  c  in the 2k-1  r o w s  of */
/*        the work array w, and  w  in turn is stored in  s c r t c h , */
/*        with  w(1,1) = scrtch(1 + nc) . */
/*            also, c  being of order  n-k, we would want  1 .le. */
/*        left-2k+ll .le. n-k  or */
/*           llmin = 2k-left  .le.  ll  .le.  n-left+k = llmax . */
	leftmk = left - *k;
	index = leftmk - j + (j - 1) * kpkm1 + nc;
/* Computing MAX */
	i__2 = 1, i__3 = *k - leftmk;
	llmin = max(i__2,i__3);
/* Computing MIN */
	i__2 = *k, i__3 = *n - leftmk;
	llmax = min(i__2,i__3);
	i__2 = llmax;
	for (ll = llmin; ll <= i__2; ++ll) {
/* L26: */
	    scrtch[ll + index] = scrtch[ll + nb];
	}
	bsplvb_(&scrtch[1], &kp1, &c__2, &xij, &left, &scrtch[nb + 1]);
/* Computing MAX */
	i__2 = 0, i__3 = leftmk - kp1;
	id = max(i__2,i__3);
/* Computing MIN */
	i__2 = 0, i__3 = leftmk - kp1;
	llmin = 1 - min(i__2,i__3);
	i__2 = kp1;
	for (ll = llmin; ll <= i__2; ++ll) {
	    ++id;
/* L27: */
	    scrtch[id + na] -= sign * scrtch[ll + nb];
	}
/* L28: */
	sign = -sign;
    }
    banfac_(&scrtch[nc + 1], &kpkm1, &nmk, &km1, &km1, iflag);
    switch (*iflag) {
	case 1:  goto L45;
	case 2:  goto L44;
    }
L44:
    //s_wsfe(&io___339);
    //e_wsfe();
    return 0;
/* ompute  scrtch (.+nd) =  d  from  scrtch (.+na) = - a . */
L45:
    i__ = *n;
L46:
    scrtch[i__ - 1 + na] += scrtch[i__ + na];
    --i__;
    if (i__ > 1) {
	goto L46;
    }
    i__1 = nmk;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L49: */
	scrtch[i__ + nd] = scrtch[i__ + na] * (tau[i__ + *k] - tau[i__]) / 
		floatk;
    }
/* ompute  scrtch (.+nd) =  x . */
    banslv_(&scrtch[nc + 1], &kpkm1, &nmk, &km1, &km1, &scrtch[nd + 1]);
/* ompute  scrtch (.+nd) = change in  xi . modify, if necessary, to */
/*  prevent new  xi  from moving more than 1/3 of the way to its */
/*  neighbors. then add to  xi  to obtain new  xi  in scrtch(.+nx). */
    delmax = 0.;
    sign = signst;
    i__1 = nmk;
    for (i__ = 1; i__ <= i__1; ++i__) {
	del = sign * scrtch[i__ + nd];
/* Computing MAX */
	d__1 = delmax, d__2 = abs(del);
	delmax = max(d__1,d__2);
	if (del > 0.) {
	    goto L51;
	}
/* Computing MAX */
	d__1 = del, d__2 = (scrtch[i__ - 1 + nx] - scrtch[i__ + nx]) / 3.;
	del = max(d__1,d__2);
	goto L52;
L51:
/* Computing MIN */
	d__1 = del, d__2 = (scrtch[i__ + 1 + nx] - scrtch[i__ + nx]) / 3.;
	del = min(d__1,d__2);
L52:
	sign = -sign;
/* L53: */
	scrtch[i__ + nx] += del;
    }
/* all it a day in case change in  xi  was small enough or too many */
/*  steps were taken. */
    if (delmax < tol) {
	goto L54;
    }
    ++newton;
    if (newton <= newtmx) {
	goto L20;
    }
    //s_wsfe(&io___342);
    //do_fio(&c__1, (char *)&newtmx, (ftnlen)sizeof(integer));
    //e_wsfe();
L54:
    i__1 = nmk;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L55: */
	t[*k + i__] = scrtch[i__ + nx];
    }
L56:
    i__1 = *k;
    for (i__ = 1; i__ <= i__1; ++i__) {
	t[i__] = tau[1];
/* L57: */
	t[*n + i__] = tau[*n];
    }
    return 0;
L999:
    *iflag = 2;
    return 0;
} /* splopt_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::subbak_(doublereal *w, integer *ipivot, integer *nrow, 
	integer *ncol, integer *last, doublereal *x)
{
    /* System generated locals */
    integer w_dim1, w_offset, i__1;

    /* Local variables */
    static integer j, k, ip, kp1;
    static doublereal sum;

/* *********************************************************************72 */

/* c SUBBAK carries out back substitution for the current block. */

/*  carries out backsubstitution for current block. */

/* parameters */
/*    w, ipivot, nrow, ncol, last  are as on return from factrb. */
/*    x(1),...,x(ncol)  contains, on input, the right side for the */
/*            equations in this block after backsubstitution has been */
/*            carried up to but not including equation ipivot(last). */
/*            means that x(j) contains the right side of equation ipi- */
/*            vot(j) as modified during elimination, j=1,...,last, while */
/*            for j .gt. last, x(j) is already a component of the solut- */
/*            ion vector. */
/*    x(1),...,x(ncol) contains, on output, the components of the solut- */
/*            ion corresponding to the present block. */

    /* Parameter adjustments */
    --ipivot;
    --x;
    w_dim1 = *nrow;
    w_offset = 1 + w_dim1;
    w -= w_offset;

    /* Function Body */
    k = *last;
    ip = ipivot[k];
    sum = 0.;
    if (k == *ncol) {
	goto L4;
    }
    kp1 = k + 1;
L2:
    i__1 = *ncol;
    for (j = kp1; j <= i__1; ++j) {
/* L3: */
	sum = w[ip + j * w_dim1] * x[j] + sum;
    }
L4:
    x[k] = (x[k] - sum) / w[ip + k * w_dim1];
    if (k == 1) {
	return 0;
    }
    kp1 = k;
    --k;
    ip = ipivot[k];
    sum = 0.;
    goto L2;
} /* subbak_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::subfor_(doublereal *w, integer *ipivot, integer *nrow, 
	integer *last, doublereal *b, doublereal *x)
{
    /* System generated locals */
    integer w_dim1, w_offset, i__1, i__2;

    /* Local variables */
    static integer j, k, ip;
    static doublereal sum;
    static integer jmax, lastp1, nrowml;

/* *********************************************************************72 */

/* c SUBFOR carries out the forward pass of substitution for the current block. */

/*  carries out the forward pass of substitution for the current block, */
/*  i.e., the action on the right side corresponding to the elimination */
/*  carried out in  f a c t r b  for this block. */
/*     at the end, x(j) contains the right side of the transformed */
/*  ipivot(j)-th equation in this block, j=1,...,nrow. then, since */
/*  for i=1,...,nrow-last, b(nrow+i) is going to be used as the right */
/*  side of equation  i  in the next block (shifted over there from */
/*  this block during factorization), it is set equal to x(last+i) here. */

/* parameters */
/*    w, ipivot, nrow, last  are as on return from factrb. */
/*    b(j)   is expected to contain, on input, the right side of j-th */
/*           equation for this block, j=1,...,nrow. */
/*    b(nrow+j)   contains, on output, the appropriately modified right */
/*           side for equation j in next block, j=1,...,nrow-last. */
/*    x(j)   contains, on output, the appropriately modified right */
/*           side of equation ipivot(j) in this block, j=1,...,last (and */
/*           even for j=last+1,...,nrow). */

/*     dimension b(nrow + nrow-last) */
    /* Parameter adjustments */
    --x;
    --ipivot;
    w_dim1 = *nrow;
    w_offset = 1 + w_dim1;
    w -= w_offset;
    --b;

    /* Function Body */
    ip = ipivot[1];
    x[1] = b[ip];
    if (*nrow == 1) {
	goto L99;
    }
    i__1 = *nrow;
    for (k = 2; k <= i__1; ++k) {
	ip = ipivot[k];
/* Computing MIN */
	i__2 = k - 1;
	jmax = (real) min(i__2,*last);
	sum = 0.;
	i__2 = jmax;
	for (j = 1; j <= i__2; ++j) {
/* L14: */
	    sum = w[ip + j * w_dim1] * x[j] + sum;
	}
/* L15: */
	x[k] = b[ip] - sum;
    }

/*     transfer modified right sides of equations ipivot(last+1),..., */
/*     ipivot(nrow) to next block. */
    nrowml = *nrow - *last;
    if (nrowml == 0) {
	goto L99;
    }
    lastp1 = *last + 1;
    i__1 = *nrow;
    for (k = lastp1; k <= i__1; ++k) {
/* L25: */
	b[nrowml + k] = x[k];
    }
L99:
    return 0;
} /* subfor_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::tautsp_(doublereal *tau, doublereal *gtau, integer *ntau,
	 doublereal *gamma, doublereal *s, doublereal *break__, doublereal *
	coef, integer *l, integer *k, integer *iflag)
{
    /* Format strings */
    static char fmt_600[] = "(\002 ntau = \002,i4,\002. should be .ge. 4 \
.\002)";
    static char fmt_610[] = "(\002 point \002,i3,\002 and the next\002,2e15.\
6,\002 are disordered\002)";

    /* System generated locals */
    integer s_dim1, s_offset, i__1;
    doublereal d__1, d__2;

    /* Local variables */
    static doublereal c__, d__;
    static integer i__;
    static doublereal z__, zt2, del, gam, zeta, temp, alpha, denom, ratio, 
	    sixth, entry__, factr2, onemg3;
    static integer ntaum1;
    static doublereal entry3, divdif, factor;
    static integer method;
    static doublereal onemzt;

    /* Fortran I/O blocks */
    static cilist io___355 = { 0, 6, 0, fmt_600, 0 };
    static cilist io___358 = { 0, 6, 0, fmt_610, 0 };


/* *********************************************************************72 */

/* c TAUTSP constructs a cubic spline interpolant to given data. */

/*  from  * a practical guide to splines *  by c. de boor */
/* onstructs cubic spline interpolant to given data */
/*         tau(i), gtau(i), i=1,...,ntau. */
/*  if  gamma .gt. 0., additional knots are introduced where needed to */
/*  make the interpolant more flexible locally. this avoids extraneous */
/*  inflection points typical of cubic spline interpolation at knots to */
/*  rapidly changing data. */

/*  parameters */
/*            input */
/*  tau      sequence of data points. must be strictly increasing. */
/*  gtau     corresponding sequence of function values. */
/*  ntau     number of data points. must be at least  4 . */
/*  gamma    indicates whether additional flexibility is desired. */
/*          = 0., no additional knots */
/*          in (0.,3.), under certain conditions on the given data at */
/*                points i-1, i, i+1, and i+2, a knot is added in the */
/*                i-th interval, i=2,...,ntau-2. see description of meth- */
/*                od below. the interpolant gets rounded with increasing */
/*                gamma. a value of  2.5  for gamma is typical. */
/*          in (3.,6.), same , except that knots might also be added in */
/*                intervals in which an inflection point would be permit- */
/*                ted.  a value of  5.5  for gamma is typical. */
/*            output */
/*  break, coef, l, k  give the pp-representation of the interpolant. */
/*          specifically, for break(i) .le. x .le. break(i+1), the */
/*        interpolant has the form */
/*  f(x) = coef(1,i) +dx(coef(2,i) +(dx/2)(coef(3,i) +(dx/3)coef(4,i))) */
/*        with  dx = x - break(i) and i=1,...,l . */
/*  iflag   = 1, ok */
/*          = 2, input was incorrect. a printout specifying the mistake */
/*            was made. */
/*            workspace */
/*  s     is required, of size (ntau,6). the individual columns of this */
/*        array contain the following quantities mentioned in the write- */
/*        up and below. */
/*     s(.,1) = dtau = tau(.+1) - tau */
/*     s(.,2) = diag = diagonal in linear system */
/*     s(.,3) = u = upper diagonal in linear system */
/*     s(.,4) = r = right side for linear system (initially) */
/*            = fsecnd = solution of linear system , namely the second */
/*                       derivatives of interpolant at  tau */
/*     s(.,5) = z = indicator of additional knots */
/*     s(.,6) = 1/hsecnd(1,x) with x = z or = 1-z. see below. */

/*  ------  m e t h o d  ------ */
/*  on the i-th interval, (tau(i), tau(i+1)), the interpolant is of the */
/*  form */
/*  (*)  f(u(x)) = a + b*u + c*h(u,z) + d*h(1-u,1-z) , */
/*  with  u = u(x) = (x - tau(i))/dtau(i). here, */
/*       z = z(i) = addg(i+1)/(addg(i) + addg(i+1)) */
/*  (= .5, in case the denominator vanishes). with */
/*       addg(j) = abs(ddg(j)), ddg(j) = dg(j+1) - dg(j), */
/*       dg(j) = divdif(j) = (gtau(j+1) - gtau(j))/dtau(j) */
/*  and */
/*       h(u,z) = alpha*u**3 + (1 - alpha)*(max(((u-zeta)/(1-zeta)),0)**3 */
/*  with */
/*       alpha(z) = (1-gamma/3)/zeta */
/*       zeta(z) = 1 - gamma*min((1 - z), 1/3) */
/*  thus, for 1/3 .le. z .le. 2/3,  f  is just a cubic polynomial on */
/*  the interval i. otherwise, it has one additional knot, at */
/*         tau(i) + zeta*dtau(i) . */
/*  as  z  approaches  1, h(.,z) has an increasingly sharp bend  near 1, */
/*  thus allowing  f  to turn rapidly near the additional knot. */
/*     in terms of f(j) = gtau(j) and */
/*       fsecnd(j) = 2.derivative of  f  at  tau(j), */
/*  the coefficients for (*) are given as */
/*       a = f(i) - d */
/*       b = (f(i+1) - f(i)) - (c - d) */
/*       c = fsecnd(i+1)*dtau(i)**2/hsecnd(1,z) */
/*       d = fsecnd(i)*dtau(i)**2/hsecnd(1,1-z) */
/*  hence can be computed once fsecnd(i),i=1,...,ntau, is fixed. */
/*   f  is automatically continuous and has a continuous second derivat- */
/*  ive (except when z = 0 or 1 for some i). we determine fscnd(.) from */
/*  the requirement that also the first derivative of  f  be contin- */
/*  uous. in addition, we require that the third derivative be continuous */
/*  across  tau(2) and across  tau(ntau-1) . this leads to a strictly */
/*  diagonally dominant tridiagonal linear system for the fsecnd(i) */
/*  which we solve by gauss elimination without pivoting. */


/*  there must be at least  4  interpolation points. */
    /* Parameter adjustments */
    s_dim1 = *ntau;
    s_offset = 1 + s_dim1;
    s -= s_offset;
    --gtau;
    --tau;
    --break__;
    coef -= 5;

    /* Function Body */
    if (*ntau >= 4) {
	goto L5;
    }
    //s_wsfe(&io___355);
    //do_fio(&c__1, (char *)&(*ntau), (ftnlen)sizeof(integer));
    //e_wsfe();
    goto L999;

/* onstruct delta tau and first and second (divided) differences of data */

L5:
    ntaum1 = *ntau - 1;
    i__1 = ntaum1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	s[i__ + s_dim1] = tau[i__ + 1] - tau[i__];
	if (s[i__ + s_dim1] > 0.) {
	    goto L6;
	}
	//s_wsfe(&io___358);
	//do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
	//do_fio(&c__1, (char *)&tau[i__], (ftnlen)sizeof(doublereal));
	//do_fio(&c__1, (char *)&tau[i__ + 1], (ftnlen)sizeof(doublereal));
	//e_wsfe();
	goto L999;
L6:
	s[i__ + 1 + (s_dim1 << 2)] = (gtau[i__ + 1] - gtau[i__]) / s[i__ + 
		s_dim1];
    }
    i__1 = ntaum1;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L7: */
	s[i__ + (s_dim1 << 2)] = s[i__ + 1 + (s_dim1 << 2)] - s[i__ + (s_dim1 
		<< 2)];
    }

/* onstruct system of equations for second derivatives at  tau. at each */
/*  interior data point, there is one continuity equation, at the first */
/*  and the last interior data point there is an additional one for a */
/*  total of  ntau  equations in  ntau  unknowns. */

    i__ = 2;
    s[(s_dim1 << 1) + 2] = s[s_dim1 + 1] / 3.;
    sixth = .16666666666666666;
    method = 2;
    gam = *gamma;
    if (gam <= 0.) {
	method = 1;
    }
    if (gam <= 3.) {
	goto L9;
    }
    method = 3;
    gam += -3.;
L9:
    onemg3 = 1. - gam / 3.;
/*                 ------ loop over i ------ */
L10:
/*          construct z(i) and zeta(i) */
    z__ = .5;
    switch (method) {
	case 1:  goto L19;
	case 2:  goto L11;
	case 3:  goto L12;
    }
L11:
    if (s[i__ + (s_dim1 << 2)] * s[i__ + 1 + (s_dim1 << 2)] < 0.) {
	goto L19;
    }
L12:
    temp = (d__1 = s[i__ + 1 + (s_dim1 << 2)], abs(d__1));
    denom = (d__1 = s[i__ + (s_dim1 << 2)], abs(d__1)) + temp;
    if (denom == 0.) {
	goto L19;
    }
    z__ = temp / denom;
    if ((d__1 = z__ - .5, abs(d__1)) <= sixth) {
	z__ = .5;
    }
L19:
    s[i__ + s_dim1 * 5] = z__;
/*   ******set up part of the i-th equation which depends on */
/*         the i-th interval */
    if ((d__1 = z__ - .5) < 0.) {
	goto L21;
    } else if (d__1 == 0) {
	goto L22;
    } else {
	goto L23;
    }
L21:
    zeta = gam * z__;
    onemzt = 1. - zeta;
/* Computing 2nd power */
    d__1 = zeta;
    zt2 = d__1 * d__1;
/* Computing MIN */
    d__1 = 1., d__2 = onemg3 / onemzt;
    alpha = min(d__1,d__2);
    factor = zeta / (alpha * (zt2 - 1.) + 1.);
    s[i__ + s_dim1 * 6] = zeta * factor / 6.;
    s[i__ + (s_dim1 << 1)] += s[i__ + s_dim1] * ((1. - alpha * onemzt) * 
	    factor / 2. - s[i__ + s_dim1 * 6]);
/*     if z = 0 and the previous z = 1, then d(i) = 0. since then */
/*     also u(i-1) = l(i+1) = 0, its value does not matter. reset */
/*     d(i) = 1 to insure nonzero pivot in elimination. */
    if (s[i__ + (s_dim1 << 1)] <= 0.) {
	s[i__ + (s_dim1 << 1)] = 1.;
    }
    s[i__ + s_dim1 * 3] = s[i__ + s_dim1] / 6.;
    goto L25;
L22:
    s[i__ + (s_dim1 << 1)] += s[i__ + s_dim1] / 3.;
    s[i__ + s_dim1 * 3] = s[i__ + s_dim1] / 6.;
    goto L25;
L23:
    onemzt = gam * (1. - z__);
    zeta = 1. - onemzt;
/* Computing MIN */
    d__1 = 1., d__2 = onemg3 / zeta;
    alpha = min(d__1,d__2);
    factor = onemzt / (1. - alpha * zeta * (onemzt + 1.));
    s[i__ + s_dim1 * 6] = onemzt * factor / 6.;
    s[i__ + (s_dim1 << 1)] += s[i__ + s_dim1] / (float)3.;
    s[i__ + s_dim1 * 3] = s[i__ + s_dim1 * 6] * s[i__ + s_dim1];
L25:
    if (i__ > 2) {
	goto L30;
    }
    s[s_dim1 * 5 + 1] = .5;
/*  ******the first two equations enforce continuity of the first and of */
/*        the third derivative across tau(2). */
    s[(s_dim1 << 1) + 1] = s[s_dim1 + 1] / 6.;
    s[s_dim1 * 3 + 1] = s[(s_dim1 << 1) + 2];
    entry3 = s[s_dim1 * 3 + 2];
    if ((d__1 = z__ - .5) < 0.) {
	goto L26;
    } else if (d__1 == 0) {
	goto L27;
    } else {
	goto L28;
    }
L26:
    factr2 = zeta * (alpha * (zt2 - 1.) + 1.) / (alpha * (zeta * zt2 - 1.) + 
	    1.);
    ratio = factr2 * s[s_dim1 + 2] / s[(s_dim1 << 1) + 1];
    s[(s_dim1 << 1) + 2] = factr2 * s[s_dim1 + 2] + s[s_dim1 + 1];
    s[s_dim1 * 3 + 2] = -factr2 * s[s_dim1 + 1];
    goto L29;
L27:
    ratio = s[s_dim1 + 2] / s[(s_dim1 << 1) + 1];
    s[(s_dim1 << 1) + 2] = s[s_dim1 + 2] + s[s_dim1 + 1];
    s[s_dim1 * 3 + 2] = -s[s_dim1 + 1];
    goto L29;
L28:
    ratio = s[s_dim1 + 2] / s[(s_dim1 << 1) + 1];
    s[(s_dim1 << 1) + 2] = s[s_dim1 + 2] + s[s_dim1 + 1];
    s[s_dim1 * 3 + 2] = -s[s_dim1 + 1] * 6. * alpha * s[s_dim1 * 6 + 2];
/*       at this point, the first two equations read */
/*              diag(1)*x1 + u(1)*x2 + entry3*x3 = r(2) */
/*       -ratio*diag(1)*x1 + diag(2)*x2 + u(2)*x3 = 0.0D+00 */
/*       eliminate first unknown from second equation */
L29:
    s[(s_dim1 << 1) + 2] = ratio * s[s_dim1 * 3 + 1] + s[(s_dim1 << 1) + 2];
    s[s_dim1 * 3 + 2] = ratio * entry3 + s[s_dim1 * 3 + 2];
    s[(s_dim1 << 2) + 1] = s[(s_dim1 << 2) + 2];
    s[(s_dim1 << 2) + 2] = ratio * s[(s_dim1 << 2) + 1];
    goto L35;
L30:
/*  ******the i-th equation enforces continuity of the first derivative */
/*        across tau(i). it has been set up in statements 35 up to 40 */
/*        and 21 up to 25 and reads now */
/*         -ratio*diag(i-1)*xi-1 + diag(i)*xi + u(i)*xi+1 = r(i) . */
/*        eliminate (i-1)st unknown from this equation */
    s[i__ + (s_dim1 << 1)] = ratio * s[i__ - 1 + s_dim1 * 3] + s[i__ + (
	    s_dim1 << 1)];
    s[i__ + (s_dim1 << 2)] = ratio * s[i__ - 1 + (s_dim1 << 2)] + s[i__ + (
	    s_dim1 << 2)];

/*  ******set up the part of the next equation which depends on the */
/*        i-th interval. */
L35:
    if ((d__1 = z__ - .5) < 0.) {
	goto L36;
    } else if (d__1 == 0) {
	goto L37;
    } else {
	goto L38;
    }
L36:
    ratio = -s[i__ + s_dim1 * 6] * s[i__ + s_dim1] / s[i__ + (s_dim1 << 1)];
    s[i__ + 1 + (s_dim1 << 1)] = s[i__ + s_dim1] / 3.;
    goto L40;
L37:
    ratio = -(s[i__ + s_dim1] / 6.) / s[i__ + (s_dim1 << 1)];
    s[i__ + 1 + (s_dim1 << 1)] = s[i__ + s_dim1] / 3.;
    goto L40;
L38:
    ratio = -(s[i__ + s_dim1] / 6.) / s[i__ + (s_dim1 << 1)];
    s[i__ + 1 + (s_dim1 << 1)] = s[i__ + s_dim1] * ((1. - zeta * alpha) * 
	    factor / (float)2. - s[i__ + s_dim1 * 6]);
/*         ------  end of i loop ------ */
L40:
    ++i__;
    if (i__ < ntaum1) {
	goto L10;
    }
    s[i__ + s_dim1 * 5] = .5;

/*        ------  last two equations  ------ */
/*  the last two equations enforce continuity of third derivative and */
/*  of first derivative across  tau(ntau-1). */
    entry__ = ratio * s[i__ - 1 + s_dim1 * 3] + s[i__ + (s_dim1 << 1)] + s[
	    i__ + s_dim1] / 3.;
    s[i__ + 1 + (s_dim1 << 1)] = s[i__ + s_dim1] / 6.;
    s[i__ + 1 + (s_dim1 << 2)] = ratio * s[i__ - 1 + (s_dim1 << 2)] + s[i__ + 
	    (s_dim1 << 2)];
    if ((d__1 = z__ - .5) < 0.) {
	goto L41;
    } else if (d__1 == 0) {
	goto L42;
    } else {
	goto L43;
    }
L41:
    ratio = s[i__ + s_dim1] * 6. * s[i__ - 1 + s_dim1 * 6] * alpha / s[i__ - 
	    1 + (s_dim1 << 1)];
    s[i__ + (s_dim1 << 1)] = ratio * s[i__ - 1 + s_dim1 * 3] + s[i__ + s_dim1]
	     + s[i__ - 1 + s_dim1];
    s[i__ + s_dim1 * 3] = -s[i__ - 1 + s_dim1];
    goto L45;
L42:
    ratio = s[i__ + s_dim1] / s[i__ - 1 + (s_dim1 << 1)];
    s[i__ + (s_dim1 << 1)] = ratio * s[i__ - 1 + s_dim1 * 3] + s[i__ + s_dim1]
	     + s[i__ - 1 + s_dim1];
    s[i__ + s_dim1 * 3] = -s[i__ - 1 + s_dim1];
    goto L45;
L43:
/* Computing 2nd power */
    d__1 = onemzt;
/* Computing 3rd power */
    d__2 = onemzt;
    factr2 = onemzt * (alpha * (d__1 * d__1 - 1.) + 1.) / (alpha * (d__2 * (
	    d__2 * d__2) - 1.) + 1.);
    ratio = factr2 * s[i__ + s_dim1] / s[i__ - 1 + (s_dim1 << 1)];
    s[i__ + (s_dim1 << 1)] = ratio * s[i__ - 1 + s_dim1 * 3] + factr2 * s[i__ 
	    - 1 + s_dim1] + s[i__ + s_dim1];
    s[i__ + s_dim1 * 3] = -factr2 * s[i__ - 1 + s_dim1];
/*     at this point, the last two equations read */
/*             diag(i)*xi + u(i)*xi+1 = r(i) */
/*      -ratio*diag(i)*xi + diag(i+1)*xi+1 = r(i+1) */
/*     eliminate xi from last equation */
L45:
    s[i__ + (s_dim1 << 2)] = ratio * s[i__ - 1 + (s_dim1 << 2)];
    ratio = -entry__ / s[i__ + (s_dim1 << 1)];
    s[i__ + 1 + (s_dim1 << 1)] = ratio * s[i__ + s_dim1 * 3] + s[i__ + 1 + (
	    s_dim1 << 1)];
    s[i__ + 1 + (s_dim1 << 2)] = ratio * s[i__ + (s_dim1 << 2)] + s[i__ + 1 + 
	    (s_dim1 << 2)];

/*        ------ back substitution ------ */

    s[*ntau + (s_dim1 << 2)] /= s[*ntau + (s_dim1 << 1)];
L50:
    s[i__ + (s_dim1 << 2)] = (s[i__ + (s_dim1 << 2)] - s[i__ + s_dim1 * 3] * 
	    s[i__ + 1 + (s_dim1 << 2)]) / s[i__ + (s_dim1 << 1)];
    --i__;
    if (i__ > 1) {
	goto L50;
    }
    s[(s_dim1 << 2) + 1] = (s[(s_dim1 << 2) + 1] - s[s_dim1 * 3 + 1] * s[(
	    s_dim1 << 2) + 2] - entry3 * s[(s_dim1 << 2) + 3]) / s[(s_dim1 << 
	    1) + 1];

/*        ------ construct polynomial pieces ------ */

    break__[1] = tau[1];
    *l = 1;
    i__1 = ntaum1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	coef[(*l << 2) + 1] = gtau[i__];
	coef[(*l << 2) + 3] = s[i__ + (s_dim1 << 2)];
	divdif = (gtau[i__ + 1] - gtau[i__]) / s[i__ + s_dim1];
	z__ = s[i__ + s_dim1 * 5];
	if ((d__1 = z__ - .5) < 0.) {
	    goto L61;
	} else if (d__1 == 0) {
	    goto L62;
	} else {
	    goto L63;
	}
L61:
	if (z__ == 0.) {
	    goto L65;
	}
	zeta = gam * z__;
	onemzt = 1. - zeta;
	c__ = s[i__ + 1 + (s_dim1 << 2)] / 6.;
	d__ = s[i__ + (s_dim1 << 2)] * s[i__ + s_dim1 * 6];
	++(*l);
	del = zeta * s[i__ + s_dim1];
	break__[*l] = tau[i__] + del;
/* Computing 2nd power */
	d__1 = zeta;
	zt2 = d__1 * d__1;
/* Computing MIN */
	d__1 = 1., d__2 = onemg3 / onemzt;
	alpha = min(d__1,d__2);
/* Computing 2nd power */
	d__1 = onemzt;
	factor = d__1 * d__1 * alpha;
/* Computing 2nd power */
	d__1 = s[i__ + s_dim1];
	coef[(*l << 2) + 1] = gtau[i__] + divdif * del + d__1 * d__1 * (d__ * 
		onemzt * (factor - 1.) + c__ * zeta * (zt2 - 1.));
	coef[(*l << 2) + 2] = divdif + s[i__ + s_dim1] * (d__ * (1. - factor *
		 3.) + c__ * (zt2 * 3. - 1.));
	coef[(*l << 2) + 3] = (d__ * alpha * onemzt + c__ * zeta) * 6.;
	coef[(*l << 2) + 4] = (c__ - d__ * alpha) * 6. / s[i__ + s_dim1];
	coef[(*l - 1 << 2) + 4] = coef[(*l << 2) + 4] - d__ * 6. * (1. - 
		alpha) / (del * zt2);
	coef[(*l - 1 << 2) + 2] = coef[(*l << 2) + 2] - del * (coef[(*l << 2) 
		+ 3] - del / 2. * coef[(*l - 1 << 2) + 4]);
	goto L68;
L62:
	coef[(*l << 2) + 2] = divdif - s[i__ + s_dim1] * (s[i__ + (s_dim1 << 
		2)] * 2. + s[i__ + 1 + (s_dim1 << 2)]) / 6.;
	coef[(*l << 2) + 4] = (s[i__ + 1 + (s_dim1 << 2)] - s[i__ + (s_dim1 <<
		 2)]) / s[i__ + s_dim1];
	goto L68;
L63:
	onemzt = gam * (1. - z__);
	if (onemzt == 0.) {
	    goto L65;
	}
	zeta = 1. - onemzt;
/* Computing MIN */
	d__1 = 1., d__2 = onemg3 / zeta;
	alpha = min(d__1,d__2);
	c__ = s[i__ + 1 + (s_dim1 << 2)] * s[i__ + s_dim1 * 6];
	d__ = s[i__ + (s_dim1 << 2)] / 6.;
	del = zeta * s[i__ + s_dim1];
	break__[*l + 1] = tau[i__] + del;
	coef[(*l << 2) + 2] = divdif - s[i__ + s_dim1] * (d__ * 2. + c__);
	coef[(*l << 2) + 4] = (c__ * alpha - d__) * 6. / s[i__ + s_dim1];
	++(*l);
/* Computing 3rd power */
	d__1 = onemzt;
	coef[(*l << 2) + 4] = coef[(*l - 1 << 2) + 4] + (1. - alpha) * 6. * 
		c__ / (s[i__ + s_dim1] * (d__1 * (d__1 * d__1)));
	coef[(*l << 2) + 3] = coef[(*l - 1 << 2) + 3] + del * coef[(*l - 1 << 
		2) + 4];
	coef[(*l << 2) + 2] = coef[(*l - 1 << 2) + 2] + del * (coef[(*l - 1 <<
		 2) + 3] + del / 2. * coef[(*l - 1 << 2) + 4]);
	coef[(*l << 2) + 1] = coef[(*l - 1 << 2) + 1] + del * (coef[(*l - 1 <<
		 2) + 2] + del / 2. * (coef[(*l - 1 << 2) + 3] + del / 3. * 
		coef[(*l - 1 << 2) + 4]));
	goto L68;
L65:
	coef[(*l << 2) + 2] = divdif;
	coef[(*l << 2) + 3] = 0.;
	coef[(*l << 2) + 4] = 0.;
L68:
	++(*l);
/* L70: */
	break__[*l] = tau[i__ + 1];
    }
    --(*l);
    *k = 4;
    *iflag = 1;
    return 0;
L999:
    *iflag = 2;
    return 0;
} /* tautsp_ */

/* Subroutine */ int SmartMathLibrary::PolynomialPack::ManagedPolynomialPackProvider::titand_(doublereal *tau, doublereal *gtau, integer *n)
{
    /* Initialized data */

    static doublereal gtitan[49] = { .644,.622,.638,.649,.652,.639,.646,.657,
	    .652,.655,.644,.663,.663,.668,.676,.676,.686,.679,.678,.683,.694,
	    .699,.71,.73,.763,.812,.907,1.044,1.336,1.881,2.169,2.075,1.598,
	    1.211,.916,.746,.672,.627,.615,.607,.606,.609,.603,.601,.603,.601,
	    .611,.601,.608 };

    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i__;

/* *********************************************************************72 */

/* c TITAND represents a temperature dependent property of titanium. */

/*  from  * a practical guide to splines *  by c. de boor */
/*  these data represent a property of titanium as a function of */
/*  temperature. they have been used extensively as an example in spline */
/*  approximation with variable knots. */

    /* Parameter adjustments */
    --gtau;
    --tau;

    /* Function Body */
    *n = 49;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	tau[i__] = (doublereal) i__ * 10. + 585.;
/* L10: */
	gtau[i__] = gtitan[i__ - 1];
    }
    return 0;
} /* titand_ */